{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Recurrent Variational Autoencoders\n",
    "Autoencoders for input reconstruction using the GConvLSTM for encoding.\n",
    "\n",
    "### Recurrent decoder model\n",
    "The first model uses a vanilla LSTM cell for decoding. <br />\n",
    "The hidden state from the encoder RNN is used to predict the last time-step of the input <br />\n",
    "and as the initial hidden state of the decoder RNN. The decoder then reconstructs the time series <br />\n",
    "in reverse using the last recostructed value as input for each step. <br />\n",
    "<br />\n",
    "The loss is chosen as the MSE of the reconstruction with the input scaled by the KL-Divergence of $P(Z;X)$: <br />\n",
    "\n",
    "$\\frac{1}{n}\\sum_{i=1}^n(X_i-\\hat{X_i})^2 + \\gamma D_{KL}(q_\\Phi(\\mathbf{z\\mid x}) \\parallel p_\\theta(\\mathbf{z}))$ <br /><br />\n",
    "\n",
    "### Linear decoder model\n",
    "The second model uses as linear probabilistic decoder parameterizing <br />\n",
    "both mean $\\theta$ and standard deviation $\\Phi$ of $P(\\hat{X};Z)$. <br />\n",
    "\n",
    "The loss is the typical VAE loss: <br />\n",
    "\n",
    "$L_{\\theta,\\Phi} = -\\log (p_\\theta(\\mathbf{x})) + D_{KL}(q_\\Phi(\\mathbf{z\\mid x})\\parallel p_\\theta(\\mathbf{z\\mid x})) = -E_{\\mathbf{z} \\sim q_\\Phi(\\mathbf{z|x})}(\\log(p_\\theta(\\mathbf{x\\mid z}))) + D_{KL}(q_\\Phi(\\mathbf{z\\mid x}) \\parallel p_\\theta(\\mathbf{z}))$\n",
    "<br />"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3d6630ae",
   "metadata": {},
   "source": [
    "### GConvLSTM encoder - LSTM decoder"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "9f33a8fd",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 200/200 [18:43<00:00,  5.62s/it]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEWCAYAAAB8LwAVAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAABiDElEQVR4nO2dd3gc1fWw37N91bstWe6922CbDqab3kIChFBCAslHSc8vlRAICQkJBAKhBEiAECAQigFTjcEGDLj3JnfJ6tKqbt/7/TGj1UpWWdmWtLbu+zx6ZubOvTNHs7v3zD3n3nNEKYVGo9FoNO2x9LcAGo1Go0lMtILQaDQaTYdoBaHRaDSaDtEKQqPRaDQdohWERqPRaDpEKwiNRqPRdIhWEBrNEYKIfCQiSkSu629ZNEcGWkFoEh4R2WV2fBf3tyyJgIjMNZ/HrnanXgYeADb2vVSaIxFbfwug0Wg6R0TsSqlgPHWVUg/1tjyagYUeQWgOe0TkEhFZJiINIrJbRB4WkQzznENE/iEiZSLiF5G9IvKGeU5E5Pdmmd+s866IZHdyn2QRuVdEtotIo4isFpFvmOdGiUhERKpFxG6WDTff9KtNOWwi8lMR2SQiTSKyUURujLn+HWb9l0XkvyLiBb7eToa5wCLzsOX6yjzXxsQkIv8yj58RkbdFxCsi75ly/c+UYamIjIy5/hQReUtEKkSk0qw37BB8TJrDEK0gNIc1InIu8Aowzdw2AP8PeMGscg3wLaAKeBJYARxvnjsd+DkQNs8tBqYCqZ3c7p/Aj836/wXGAs+IyJVKqR3Ap0AWcIZZ/6vm9kWlVAC4C/gjIMB/ABfwmIhc2+4+lwGjgWeBsnbnioH/mfsNGCalBzqRt4WrgUagBjgTWANkADuAY025EJHB5jM4E/gE+Ai4FHhXRJzd3ENzBKIVhOZw5xZz+3ul1LXAXCAEnC0i4wC7eX4d8BxwPZBnlrWcK8Lo8G8BhgB72t9ERPKAy83DM5VS3wR+YR7fam6fMbdfM7ctCuIZEZEYWT8DmoD15vF3291uB3CMUupGpdQ7sSeUUkVAiympRin1faXU99vL244PlVKXA/8wj70YSqBF/pnm9htAJsbz2AOUAJXABODUbu6hOQLRPgjN4c4Ic7sJQClVJSJVwGBgOEanPRe4CLgCUMAHInIJ8B7wd4yOscVssxy4ECjt5D5epdRuc3+zuR1ubv8LPAhcLCITgFnAVqXU5yKSC6SY9a5vd+0x7Y6/VEqFuvvHe8Amc+sxt0VKqYiINJjHyeZ2hLmdaP51JaNmAKBHEJrDnV3mdgKA6T/IMct2AyGl1NeANIxO7wOMt+dLASvGW30GRgf4DEan/q0u7uOOscmPj7kPSqk64HUgHXjcPNcyqqjCGDUATFdKiVJKMH6Ds9rdy9/N/xw2t/H+fsPdHLewy9y+2iKfKWM+hglOM8DQCkJzOPFHEfk85u9k4GHz3C9E5F8YdnMb8L5SaitwpYhswvAffA/DxwDG2/TxwE4M09MPgRNizrVBKVWBMY0U4H0ReQr4vXkcO3uoRSGchDFaedZsr2Jkfc90nD+PYU66o0dPAfaa20IReUJE/q+H7TvjOYz//RLTWf+YiHxg3m/QIbqH5jBCm5g0hxPj2h1nKaVeE5GvAj8DvoLhiH0Mw/kMsAXj7f1cDOdzKfA74E2MUcM2DGd1hlnvUVrf/tvzTYzO8hIMP8N24D6l1H9i6ryL4VgeDHyslIr1Z/wKqAauw3Ac1wOrgBfj/P8BUErtEpE/Y4x0bgA2YDi/Dwql1D4ROQW4G5gDnIjhi3gY49loBhiiEwZpNBqNpiO0iUmj0Wg0HaIVhEaj0Wg6RCsIjUaj0XSIVhAajUaj6ZAjahbTvHnz1DvvvNN9RZOamhoAsrKyekukAyZRZdNy9QwtV89JVNmOYLmksxNH1AiiqkrPxNNoNJpDxRGlIDQajUZz6NAKQqPRaDQdohWERqPRaDrkiHJSazSaviEYDFJcXIzP5zvk145EIgCUl5cf8msfDIe7XC6Xi8LCQux2e5f1YtEKQqPR9Jji4mJSU1MZMWIERqqLQ0coZEQ6t9kSq3s6nOVSSlFdXU1xcTEjR47stF57tIlJo9H0GJ/PR3Z29iFXDpreQUTIzs7u8YhPKwiNRnNAaOVweHEgn5dWEBqNRnMYU7/bi6eomZCvszxQB06vKQgReUpEKkRkfSfnfyIiq82/9SISFpEs89wuEVlnnlveWzJqNJrDF6vVyowZM5gyZQoXXHABHo+n32T56KOP+Oyzzw7Z9V577TU2btwYPb799tv54IMPOqwbbAoTbAzRG5kbenME8S9gXmcnlVL3KqVmKKVmYCR3+VgpVRNT5VTzfPt0jBqNRoPb7Wb16tWsX7+erKwsHn744e4b9RJdKYgWJ3JPaK8g7rzzTs4444wO60ZCxiwmi7XHt+mWXnPHK6UWi8iIOKtfCTx/sPcMhULRuCTx0J9vHN2RqLJpuXrGkSpXJBI5oI4vHsLh+E0lLTLMmTOHdevWEQqF2L59O7fddhuVlZUkJSXx6KOPMmHCBMrLy7n55pvZsWMHAA899BDHH388999/P08//TQA119/Pd/73vfYtWsXF1xwAccffzyff/45BQUFvPTSS7jdbv72t7/xj3/8A6vVyqRJk7j77rt59NFHsVqtPPvsszzwwAM89dRTuFwuVq9ezfHHH09aWhopKSn88Ic/BGDGjBm89tprjBgxgmeffZb77rsPEWHq1KncdNNNzJ8/n48//pi77rqL//73v9x9992cd955XHbZZXz44Yf89Kc/JRQKMWvWLO648Q84nU5Gjx3NN675Bm+++SbBYJAXXniBCRMmtHlekUhkvz6yqxhO/T5fS0SSMEYat8QUK4y8vQp4TCnVWQpIRORG4EaAwsLC3hRVo9F0wAT7/3rlupuDl8VVLxwOs2jRIq6//noAvvvd7/Lwww8zduxYvvjiC2699Vbef/99fvCDH3DSSSfx8ssvEw6HaWxsZMWKFTz99NN8+umnKKU44YQTOPnkk8nMzGTbtm08++yzPPbYY1x55ZW8+uqrXHXVVdx7771s27YNp9OJx+MhIyODG2+8sY0CeOqppygpKWHJkiVYrVbuvPPODmXfsGEDf/jDH1i8eDE5OTnU1NSQlZXF+eefH1UIsfh8Pm644Qbeffddxo0bx7XXXMczL/+TG7/xHRDIzs5m2bJlPPLII9x33308/ninXWdc9LuCAC4APm1nXjpRKVUiInkYCeI3K6UWd9TYVB6PA8yaNUsdSETDRIvOGEuiyqbl6hlHmlzl5eW9vh6gu+t7vV5mzZpFSUkJEydOZN68eXi9XpYuXcqVV14Zref3+7HZbCxatIhnn30Wm82GzWbD6XTy73//m0svvZT09HQALr30UpYuXcqFF17IyJEjmTXLsHDPmjWLvXv3YrVamTZtGtdeey0XX3wxF198MTabDYvFgsViicpssVj46le/itPpjB7Hnm/5/z5Y8CHnzr2QnMxcbDYreXl50fpWq7XN9axWK9u3b2fkyJFMmjQJgG9cdQ0P/uVBbrruuwBcfvnl2Gw25syZw+uvv77fM7RYLD36zBNBQVxBO/OSUqrE3FaIyKsYCdQ7VBAajaZ/KVJfO6TXi9d01eKDaG5u5uyzz+bhhx/muuuuIyMjg9WrVx+0HC2dOxgO8Ra53nrrLRYvXswbb7zB3Xffzbp16zpsn5ycHN232WzRFc9AdD1CsDmECkXwe4LYBrV1IoQDEZRSXU5PVeY1xSJtZI6V92Do12muIpIOnAK8HlOWLCKpLfvAWUCHM6E0Go0mKSmJBx98kL/85S8kJSUxcuRIXnrpJcBYQbxmzRoATj/9dB555BHAMEvV1dVx0kkn8dprr9Hc3ExTUxOvvvoqJ510Uqf3ikQi7N27l1NPPZU//vGP1NXV0djYSGpqKg0NDZ22GzFiBCtXrgRg5cqV7Ny5E4DjjzqRNxfOp6K8EmjN7eCyutm3uZyabY2EA60+mfHjx7Nr1y6KiooAeO4/z3HcUScg1t5Zk9Kb01yfB5YC40WkWERuEJHviMh3YqpdArynlGqKKRsEfCIia4AvgbeUUvFnAdJoNAOOmTNnMm3aNJ5//nmee+45nnzySaZPn87kyZN5/XXj/fOBBx5g0aJFTJ06laOPPpqNGzdy1FFHcd111zFnzhyOOeYYvvWtbzFz5sxO7xMOh7n66quZOnUqM2fO5LbbbiMjI4MLLriAV199lRkzZrBkyZL92l122WXU1NQwefJkHnroIcaNG4dSirHDxnPb9T/g3K+czfTp06M+jIvOvJRH/v0Qp1xwAl8uWI3PEyDkD+NyufjnP//J5ZdfztSpUxERvnHZdUgv9eSiemPybD8xa9YstXx5/MsmEjVDFCSubFqunnGkyrVp0yYmTpx4KEWKcjjHPOoJ4WCEijUeAOwpNnImpEXPla/1EAlEsLqshP1hUODOcZAxIqXNNer3NtNU7sORbiNzVAoWa9eaopPPbWBklNNoNJr+IOQL49nVRKApfrt/2N/qk4gEY/ZDESKBCAikDnGTMtgFQLB5/+m/LWsgeqsn1wpCo9FoDoKQL0z1lga8VX5qtjUQ8se3jiPWtxAJtVpygl6j3Gq3IAJWp+G8DvsMp3UsLe0OOx+ERqPRDARqtzcaIwABFVLUFjXu15F3RDjQOmpQYYWKGG1C5kjB4jC6Z4tNEAuoiCISbKcgzJGHRY8gNBqNJrEIByOEvGEQSBuahFiFkDfcoTlov7YxCiL2OHYE0YLFYYwigs2GCaupwkf93mbCLSMIix5BaDQazSElElGEg5HuK3ZCqKUzd1iwOizY3GZH3ti9LyLWBwFETVMtSsDqau2ereZoItgcJhKORJ3TEVOpaBOTRqPRHGJqtzZQscbTxh/QE1pCbFvMt32b0+zI43BWt4wYxCbmtSI0lDRHlY7N2bpwLnpdb5hAQ8gIRtSC0GvTXLWC0Gg0hyUpKa1TPhcsWMC4cePYvXs3d9xxB3/+85/3qx8JK7y1ASLhCHfccQdDhgzh5AtP4IRLZ3PppZe1iZ76rW99q81xZ7R05i0KosWh3GIm6ooWBWFzGW0a9jbTWOoDBVanBYutdVTQ4o8IecP464NtriMWgV5K3pRYE401Go2mhyxcuJDbbruNd999l+HDh3dYx1cXoG53M5FABGe6HYBbbrqV6869EYD3Vr7Faaedxrp168jNzeWJJ56I697Rt31Hi4IwO3JfGKUUkUgEq3X/ONyRUAQVViDG6CDY2HouKdeJI7Vt19xiYgr7w/jr+m7tmh5BaDSaw5bFixfz7W9/mzfffJPRo0d3WCfkD1Nb1Bi11/vrgoQDkTYZ2C49/zLOOuss/vOf/wAwd+5cli9fzqOPPspPfvITwJgx9K9/Ps0ttxiBp5999lnO+MqpnHHVKXz/l7cRDoexWIUxJw/jt/f9munTprN06VKefPJJxo0bx5w5c/j2t7/NLbfcQnOln+raKr79f9dx2sWncM41p/Plmi+wua3c+8g9fOe2mzjr/DOZOGM8Dz/2EBarYHVaeOnNF5h76QmccdXJ3Hq7EZSiqqqSq795FcccewyzZ8/m008/PWTPV48gNBrNQXG3vNgr1/1lN0EA/X4/F198MR999NF+eQ9aUEpRt7sZFNiSrFjtFvx1QfwNQSz21u4vElQcddRRbN68uU37yy67jOOOPY5ffOc3+GoDPP/M89zx+9+wadMmXnzhRV5/YgF2u53f/P3nvPDS83z9iqtp9jYxc8rR3PWL31NRW85dd13N8uUrSHYlc9Y5ZzJp7GQaSrz8+i+/4OYbb+HEk05i84ptXHXrV1j+0SoAtmzdwrtvvEdDYwPTZ0/lxm/exB7PDh546j5ef/JtBhXkUd/kAeD2+3/Bzd+9hTPOPY3ikmLOPvtsNm3adBBPvhWtIDQazWGJ3W7n+OOP58knn+SBBx7osI6/LkigPohYICnHiVgEf32QsD+CirGfREKKQFMQX12wdXUykJuby9D8YXy6+FNGDh3Nth3bOHb2cTz2xKOsXLmSc645A7FAIOwnL9cI1W21WjnvtAsIecN88fkXHD/rBKiw4fE2c+nFl7F++QYAPln2MTtKtqGU4Y9obGogIEaU13POOgen04nT6SQ3N5fyinIWL/2YSy+5jEEFebjS7aTk5+OvD7Lky48p2rMVy6+Nf6i+vp7GxsY2PpoDRSsIjUZzUHT3pt9T4g1TbbFY+O9//8vpp5/O73//e37xi1+0Oa+UMpy+wL1P/IEPlrwHwJK3lxrTQsWw7YcDESKhCCu+XMmUsdNoqvRHrxFoCnHBaZfwxgevM3bMOM6Zex6B+hBKKa66/Ov86JqfY0+xRcNhALicLuwue9TJHPK1rosINASJmL4HheLj95fgcrmMWU8i2JMMf4XD6Yhez2qxEg4bz8RiFVIL3K33SrcTUYqFCz4ie0hGt7GYeor2QWg0msOWpKQk3nrrrWgE11gCjSGCTSHEItx11+/4YskyvliyDEeKDVe6HXe2A3eOkT/h9bdfY9GnH3LJ2ZfhqwlEr9FU5uOcuefx3pK3ef3dV7jorEsINIQ47bTTeOX1V6iqqcTmtFBTW8PuPbuNRgLpw5NIHeJm9tGzWLryMzz1HkKhEK+98RpgTGE9/dQz+PvjRh5te7KNjdu7zmow96RTeeX1/1FdUw1ATa0RcPH0U0/nsScfidY7FLkwWtAjCI1G0++oiDrg1cBuSWH+y29wxrmnk5ubC8Dvfvc77vvz/aAUiLBj88792v3tkb/x/IvP01jfyIQxE3np76+SnZlDyBtGRRSBphA+AmSkZzBh/EQ2bd3EzMlHE2gKMWbkOH560y+48tavgA0cdjv33/sAw4e1nUU1avJIvnfDDznvm2eSkZLBmBFjSUtJw+qy8Jc/3sf3f/I9Zp9wNKFwiBOPO5G/3f9wp//npImT+L8f/YyzzjsDq9XK9Gkz+Mffn+DeP/yZH/74B8yYOYNQKMTJJ5/Mo48+ekDPsj063DeJF4oZElc2LVfPOFLlOlThvuuLm2mu8qNCirRhSSTnuXoUVrslZLbYhEHTM6LZ11REUbaqFhSkDUuKThNtTySkqNvVtF+5O9uBvyFEJBDBkWYnOc9JOBymfpcXlBGeO9gYwpluJynX2cGVW2nxB9TsrOfa267myguv4itXfgVHyiEKGx4JEwkpkrNdhzzctx5BaDSafkFFFE3lvuiqYG91gOQ8V9eN2hE2w1OokCISUljtRl8XDkRAGauUO1MO0HmICm+1YWayOi0k5Zj+ADGOw76IEUpDiK6p6Irf/fEuFn30Ic3NXk6ePZd5p54X9TUkOlpBaDSafiHkC7cJGdE+eF08xMYzCnnD0QB30RAYtq7fqMUMU6HMy9jcVqxOCyGfUZCU62hj+rKnWAn7IlgcFtyZ9i6VTwv33PXHqHwNJV6sLmuvBdc71GgFodFoDgilVNSkcyC0zOyxua2EvGEiwUg05HW8hGIVhC+MM81u7hvlLSOKrhCrBRUxw2bbhKSczk1GjjQbjmQ7Fpulx9EtbG4rKQXuNiE0+pIDcSfoWUwajabHuFwuqqurD6jTaSEatTQm7lDs6uZ4iB11xLZtMT11N4KAtmameKKitiTyORDsSda4Rh2HGqUU1dXVuFw9M+H12ghCRJ4CzgcqlFJTOjg/F3gdaJle8IpS6k7z3DzgAcAKPKGUuqe35NRoND2nsLCQ4uJiKisrD/gajWU+wv4wjoCdkC9MJBChcqMTq9tMlBNHFpzGch9hUzHYGq0kN7nalDsCdmyerq/TEnoDwOGzYavr2D/QogwPZtTUGxgxnxTOSnuXpiuXy0VhYWGPrt2bJqZ/AQ8Bz3RRZ4lS6vzYAhGxAg8DZwLFwDIRma+U6j60okaj6RPsdjsjR4484PYqovjznFcINIY457Gj2fH6PrYvKOXE2ycx9Xv5QHwzrB46983oLKTMMSn8v23nAfDgmfNpKPFyxv0zSBnd9YridQ9sZdfCCgCO+/kECo/r+L4NvjoAUl3p8f2TfUSDr47mCj+Trx2LI6V7p3lP6DUFoZRaLCIjDqDpHKBIKbUDQEReAC4CulUQoVAoOn0vHjwezwGI1zckqmxarp5xuMkVqA+x8v6dTL6hkNRCd4d1umLPh1WkFLjImtB1p+zZ3kSgMYQj3UokM4DdWL5A6foqhnqS47pXOGgkzmmhvriZ6upqQs0RGkq8iBUkLxjt2DtDklvNVJIV6rR+U6Cxw/L+pinQiDfsp6a2Bkeg5wqiK0Xc3z6I40RkjYi8LSKTzbIhwN6YOsVmWYeIyI0islxElldXV/emrBrNEc+Gf+5l1YO7+OLuoh63rd3WxFtfXcWCK1d165uoWtsAQEqhEwGS8o2ppLWbm9j4WCl1W73d3q+pxIcKKxzpNmzJxuyi5jI/dTsNpeHKcbRJ29kZ9pjQ2u5cRxc1Bx79OYtpJTBcKdUoIucCrwFje3oRpdTjwONgLJQ7kIU/ibaIKZZElU3L1TMOF7kad2wFoHpNU/Rc8WdV1GxrYNq1XZuUdiyvBgUNe3xYPA4yR6d2Wrd+0x4AMoelGSab4Q5gNzUbm6jZ2MSul6uZvnUiKgLLHtxK2YpaTv/LDFLzW0c1dauMxDnJuS7EAjVbGwkWC+F6QymkDnLHZQ5KyzKUkT3FFtfnlGgmJgCx2snKzDp8TEzdoZSqj9lfICJ/F5EcoAQYGlO10CzTaDS9TNUm42fp2dFIJBTBWxPg+XkfE2gIMfioTPKmZkTrbn+3lD0fV3Lyb6dgtVvY/VFF9FzRglJm39q5gij53BjtZ4036iQPaju1tG6Hl49+sY6iBaVUrjdMPvuW13DNktNJGWQ4oj07DZOPO9uBzWWlZmsjFRvqCJnTZ5PiXHTXstjNFceit4FGv5mYRGSwmNMBRGSOKUs1sAwYKyIjRcQBXAHM7y85NZqBgoooqjcbpp9IUFG5oZ7Fv1lv5EAGipcanbpSisV3rOeFeYv57A+b2PzyXpRS7PmodUbTrg8r2ly7fK2H8jW1AIQDYcpWGPt5U423cYvNQtrwJAAyxhvbpX/aTOX6OpLynCQPdlG7rZFXLm9NhuPZYTin3dlO0oYabSrWeihf4wEgLU4fStbYVNw5DvKmJ97IoL/pzWmuzwNzgRwRKQZ+A9gBlFKPAl8BvisiIcALXKEMw2VIRG4B3sWY5vqUUmpDb8mp0WgM6oubjbDTJuue3cWqx3dEjyvWegDD5LTkt60/ye3vlpE3LYPmKj9iE1RIUfxZVfR8yB/mmRMXEmgIMe6iIcy8cRQhX5jkQc5oNFWAY388Hs+OJtJnO1j6syIa9/jIn5PF0TePRkR4+6YV7F1SRdmqWlyZDlY9th2AtKFJ0XDb5as8BJuMEUTm2M5HMLG4Mh2c/9ScHj6tgUFvzmK6spvzD2FMg+3o3AJgQW/IpdFoOqbFvNTCsge2oiIKV6YDX22Ayg2Gqaf4M2MkkZLvorHUR8lnVVHz0uAZGVRtbqC5wk/tjkYyR6Xg2dEUHYVsfb2EnR+UA5A+MqXNmoL04cmkD0+mwVfHnN+MJLzPSvbkNKzmYrcRZ+RR9EYpH/1qHY37vHhrAuRMTmPEGXnR3NDVmxsIByJYbELWmPhmQ2k6p79nMWk0mgSh2lQQLTb5SMhIbDPl6mGA4QQGKFthTCUfdkouYhVqihpZ/5yRCyFrfCq5k9MA2P52qdFum2G2Sh7sxOqyREcpGSM778AdqTbypmdElQPAmHMLjOsuKKV8tYekPCezvzcOm9OKK8OBO9sRXfCWku/G5taRhA4WrSA0mgFIY4mPz+/ahq+uNTlOywhi8FGZ0bK8aekMPTkXLNBY6iXkC1O20vAf5E5JJ2tsCigoWVqN1WVh+Nw8cqcYtvxdHxojhRYFkT0hjUlfGxa9ds6k+ExALaQOcTP4aEO2zDEpnPjrSW0yuWWOaV17kTa052s4NPujVaxGMwBZ9eAu1j+xF5fLzezbxvLZHzZFVxMPPiqDfV9WE2wKM+qswdicVlLz3TSUeNn7SSU12xpNE04KORPToo7tYSflkjzIFVUQxZ+amc+2GSOP5DwXYy8sYO8nlQSbwmRPSOux3Mf+eDxlq2oZfFQm9qS23Vfm6BT2fWGMblJNp7Xm4NAKQqMZgNRsNjrtPR9VEGoOsfxv26LnMkYlc8wPx+PZ1cSQ47IBI+lOQ4mXtc/sAiClwI0tyUbO5DS2vFqCxW5h/CXGetaMEcnYk600lfvw7GykZquhQFKHuLDaLZzx5xkoUXHFWmqPPdnG0BNzOzyXGRNSI2tM1yu5NfGhFYRGcwTgrQ2w8pEidrxbRsGcLE76zeQuF03VbTdWG1es9eCvNxacpQ11kzo0iZQCN2lDk8mf3bpoLG1YEiVLq9nyakn0GGDwzEyGzc0lrTCJ1EKjTKxCzuR0Sr+sYfs7ZVETU+rQ5Oh56TyJ2QETNTEJZI7vmflK0zFaQWg0RwCvfu0zdr5v2Pz3LK5k8/+KuX7ZmSRlO2ko9bLqse0MOTab0fPyCTQGaSr1A0ZOhqqN9VgcFk7703TsyR13CZmjjM492Gg4mNNNBWGxWzjmh+P3q583xVAQRW/to6HYi1il1/0CrkwHU68dTrAxjDtDh8w4FGgFodEc5jSWedn5QTkWuzDl68PZ/k4Znp1NrHtmFxmjUnj9qqUEm8M4M+z8sOpiqk2TTyyZo5I7VQ4ABcdkM+b8fLYvKEUpojOVOiPXXABX9JYxk8md7cDm7P00mxMuG9p9JU3cdKkgzNDb5wFblFJb+kYkjUYTD4t+uZY9H1cy/uIhoCB7QirjLy3E5ray8pHt7Hi3lIZiL8HmMGIBvydI6fJaPDv2j0qaPb7rDl8swswbRzPuwiE0lHrJmdz1quOMkcmkDU+ifrdhykrO6zxLmyZx6dJLpJQKA08Cx/WNOBqNJl4++/0mij+tYuFP1gAwaFqGsZ1ubHcvqqRyQz32JKsxVRXYOr8kOoJIG9Nq8smNM8xE8mAXg2dmdptTWSzCrFtaY2+2n3GkOTyIZxrBc8B1IjJZRLJa/npbMI1G0zmx6xdaGHJ8DgDJ+S6Scp3RRWO5U9MpmGP8ZHd+UE7NFkNB5M1Kw5lux5luJ2/KoY9DlD0+lXEXG4vbBs3MOOTX1/Q+8aj12wAFrI0pU3G21Wg0B4FSiufP/phAY4hrlpyGxWq803l2NrWplz48KZrgR0TIm5beuq7h6EzypmWAQPmqWoJm2Iu0EU7OuH8GkUAEm6t3/APTrh/JmPMKcOdop/HhSDyd/GIMhaDRaPqAcDDCc6ctwp5kY+4fpkZnJ5Wt8lAwyxgJeLYbfoTM0SlkjEk2zD4xcY0Gzchg18IKrE4Lhcdn40yzkzEqGc/2puiK6dThbpJyetc3ICIkD4ov7LYm8ehWQSil5vaBHBqNxmTLq8Xs/cSIhhr0tkZX3bWwPKogas1Q12nDkph18/55tvJnZ5E9IZWM0Sk404y39+Fz8/Bs3wmAI8OGK1u/1Wu6plsFISLpwAPAOWbRW8APlFJdJ3rVaDQHROyq5r1LqmL2K+H/JgJQa44g3J2MAOxJNk770/Q2ZeMuGkL+7CzKV9ViyQv3wlI1zZFGPE7qB4FrgID5dx3w194TSaM5ctjywj7e//Zags2h7isDZatr2ftJFVaXhZYeXMxf6b5lNdF6LVNVU4f0zHyTWuBmzHkF5E7teRwkzcAjHgVxDvAnpdRQpdRQ4F6MtREajaYbVj24i6JXytn00t646q9/dhcAQ47Njs48GnRUJrYkK80Vfur3GqallmxqqUN0UDpN73Eg4b61w1qjiZOmfUZIi21v7tvvnFKKrfNLqFjviZbt+9IYJQyakcnUa0Yw6KhMJlw2lGwzttDOD8qJhCLU7TYVRJxpNTWaAyEeBbEA+ImI7BGRPcBPMPwQGo2mCwKNwWgmtd0fVVC+1sOXZpY2pRTv3LyCly76hFe++hlg5IQuW+UBIGdiKmlDkzj5jsnkTk4jZ6JhEtr5QTn1e5uJhBTODDuOLsJjaDQHSzzfru9jKJIWJ/WzwA96SyCN5kihocQb3fdWBXj2pIX460O4sx3U72lm5SNGTuXqTQ0opaje2kCwKYQr007y4La+hfzZWWz4zx62v13KtGtHApCUq8NXaHqXeGIx/Rr4p1Lqmp5cWESeAs4HKpRSUzo4/3Xg/zBccQ3Ad5VSa8xzu8yyMBBSSs3qyb01mkQgVkEA+OuN0cSap3ZSub7tJMCmcl80lWfasKQ2axrAyNGQPMhFU7mPxb9ZBxgJeDSa3iSeWEwXA6MP4Nr/AuZ1cX4ncIpSaipwF/B4u/OnKqVmaOWgSSQ+/cNGvvxrfHErWxSExSHm1vi57V5UQXOln5QCF+kjDCdzTVFjNJVnWgfZ0ESEwhONUBoln9eAwIgz8g7un9FouiEeE9NHwO0i4gRKWwqVUq901UgptVhERnRx/rOYw8+Bwjhk6ZJQKERNTU33FU08Hs/B3rLXSFTZBrJc3qoAH/1iHWIVRl+d023AuvJtRsrN7DluHA4nGWPdlHzkwbPFiHCaMyeFpj1+2AXF68rZ+0UlAK7hFhp8+y8zyprtgv8Z+3mzU0mebO2wXjw0BfaP6JooJKpsiSyXN+ynprYGR6DzJFGdkZXVeWi9eBTE9eb2QXMrGDOZDmXwlhuAt2OOFfCeiCjgMaVU+9FFFBG5EbgRoLDwoHWMRtMp1RuMIHcqrPDVBnF3sxK5JSmPM8fGlKuM72bYr/BsaUYsUHhqFrvfNBbCNez1UrnGCIGRMTa5w+uljXSROsJFc1mA0Zfn6YVuml4nHgXx294UQEROxVAQJ8YUn6iUKhGRPOB9EdmslFrcUXtTeTwOMGvWLNWVNuyMA2nTVySqbANRrm07K6P7dp+LrKyM6LG/IUhjqY+MkclY7YYpKVhtRFNNzUkm1WVESx13ehJ736kla2wK+SPz8Azys5caqpY1EWw0kvoMHpHb6ejk9HtmEPKFcWcdGgd1i1yJSKLKlohyidVOVmZWl2lmD4R4nNRpwJtKqUWH9M7G9acBTwDnKKWqW8qVUiXmtkJEXgXmYAQN1Gj6jYq1nuh+U7mP0kAN4UCEIcdk88wJC6lYV4fVaeH8J2cz5esjoj4IV07rj9aZbufcx1vdai0zkfYuMZRP2tCkLk1X9iSbzq2g6TN600ndJSIyDHgF+IZSamtMebKIpLbsA2cB6w/1/TWanlIeoyAay3w8Net9nj5+IWuf3knFOsMXEPZHWPmoMX21RUEk5XZuimoJgx32G6ONtGF6ZbQmceg1J7WIPA/MBXJEpBj4DWA32z4K3A5kA383p/S1TGcdBLxqltmA/yil3unRf6XRHCD7ltew8YU9nPxbY2Z25ToPQ47NIRyMULWhPlqvtqjVYfn2d1cARojt8tUeqrc2EAlHaCzzAeDO61xBtA+3nTGyY/+DRtMf9JqTWil1ZTfnvwV8q4PyHcD0/VtoNL3P4tvXs/3tUvKmZVC9qZ7P7tnExf85lrxpGdEMbQCVG1uVRcvb/8SvDqVifR3NFX48O5pQYYU9xYrN2flPpX001uyJOoieJnGIR0HciY6/pBkgVKzzAFC9pZ6y1ca6hPXP7WZyu19A9eb6NscZo5LJmZRGaoGb+j3NFL1tDLYdGV3/xKx2C84MO35PEGeajbQhOraSJnGIJ2HQHX0gh0bT7/jqAjQUG36D+r3NNOw11ivs/aQqavpxpNkI1Ieo3WZMeU0pcGFzWRl/yRDEIqQNTaJ+TzPb3jCC8zm7URBgmJn8nqCxgrqbtRUaTV/SqZNaRFaKyJmm0/gpEZlgll8iIvGvRtNoDhOqN7WOChqKm6k3lYW/Lsiqx3cAMPQEYzVzyGeYlbInpHHmX2cy7BRjVXPaUGMEsOsDI01o+ujuRwQtfoi0odr/oEksuprFNAPIBFwYSYIKzHIHkHgTgTWag6Qyxgnt2dGEvy4YPQ4HIqSPSGbI8Tlt2rjS2847bz8LafBx3f9UCk/IxpVhp/CE7AMRW6PpNbob/2rfg2bAULmhNWxF3Z7m/c4PPzUXV0ZbheDKajtDKTaOUmqhm7SR3U9bHXZKXnQEotEkEt0piG8CZ2IoiltE5GJg/wzpGs0RQOw01pZXo5QCF97qADaXlZFnDiIcaPvO1D7cRkqBG7GAisDgmZmIdiloDmO6UxBnx+xfHLOvRxaaI47oCKJlIjfGiODEX00CwJFiJxKMtGnTfpqq1W4hY3QKnu2N2mSkOezpSkGc2mdSaDT9jM8ToKHEi8UupBYmUbfTSOnpynCQWthqJrLYLdiSrISawwAk5+y/CO74n02kscxH9sQ0Gv31+53XaA4XOlUQSqmP+1IQjaY/2fupEVU1Jd9Ncq6zVUFk7h/8zJlqI9QcxmITnJn7K4ikXKfO9qY5IognJ7VGc8Sz7uldAAyamdHGbNRRR+9Is0e3Fqv+CWmOXPS3WzPg8dYG2Pp6iZGl7bQ83DEzk5I6SOvpNBWEM/3QhlbWaBINrSA0A56NL+4hHIiQPT6V9BHJ0QirACn5+ysIR6phmW2/BkKjOdLo1AchIrd31VApdeehF0ej6VuUUqz+h7FKeshx2YgI7mzDrCRW6TBrXHQEkabzMmiObLr6ht8Rs68gmuGwZYqrVhCaw57iT6soW1mLI9XGiDMGAZAyyBg1uLMcHfoYWhbDpQzRuRs0RzZdKYjLze2pwCnA/Rgmqe8BS3tZLo2mT/jyASNXVeEJOThTjZFB8mAXc344Dkdyx2G6R545iKyxKaQW6sirmiObrqa5/g9ARP4M3K2Ueso8FuCnfSOeRtN71BQ1sOXVEsQqjDkvv8254XM7D30hFiFjVEpvi6fR9DvxGFGdwG9EpBBjBHE92rmtOcyJhCO8cd2XqLCi4Jgs0ofrSKoaTXviURA/Bp7ASBEK4MOI0aTRHLZ88ZctFH9ahTPDzvRvjuxvcTSahKTbkYBS6j/ACIxYTBcDI5RSL8RzcTOPRIWIrO/kvIjIgyJSJCJrReSomHPXisg28+/aeO6n0cRDJBTh8z9vAWDK14eTkq99CRpNR8RrKpoNnAYUAWeJSLw5o/8FzOvi/DkY0WHHAjcCjwCISBbwG+AYYA6GiSszzntqNF2yc2E5zZV+kgc5GXG6DrOt0XRGtwpCRL4PvAHcCgwGLgXujefiSqnFQFfZ5y4CnlEGnwMZIpKPEUX2faVUjVKqFnifrhWNRhM3G5/fA8DgWVlYbNqdptF0Rjw+iO8DLwFfMY8/4NCtgRgC7I05LjbLOivvklAoRE1N/NlQPR5P3HX7mkSV7XCXK+QNs+l/xlcr9wQ3Db66blocHE2Bxl69/oGSqHJB4sqWyHJ5w35qamtwBHq+uj8rK6vTc/G8PmUCa2KOk4COJ4j3AyJyo4gsF5Hl1dXV/S2OJsHZ/X4VwcYwKUOdpI/RM5c0mq6IZwTxJfBdc//HwInAp4fo/iXA0JjjQrOsBJjbrvyjji6glHoceBxg1qxZqitt2BkH0qavSFTZDle5Fr25CYAhc3JJc/VdavXUPrxXT0hUuSBxZUtEucRqJyszC0fKoY0PFs8I4lbAixFqYx5QimF2OhTMB64xZzMdC9QppUqBdzGc4Zmmc/oss0yjOWD89UG2vbkPBIafktPf4mg0CU+XIwgRsQLjMBzTLbkWtyilwvFcXESexxgJ5IhIMcbMJDuAUupRYAFwLsbsqGaMRXgopWpE5C5gmXmpO5VS8TsXNJoO2PJaMWF/hMyxKaTphXEaTbd0qSCUUmEReRL4iVLqXz29uFLqym7OK+DmTs49BTzV03tqNB2hlIomBcqfnYkRMUaj0XRFPD6I54DrRGQZhnkJMN7ye00qjeYQs+WVYnZ9WIHNZWXEqYP6WxyN5rAgHgVxG0aI77UxZSrOtocFfl+YTz8oZ/umesZOTmP1FzVsWl3L//1pOqPGp/W3eJqDxOcJ8O6tKwEYe1EByYP2TwKk0Wj2J55OfjGtOSCOOJZ+WMUvrl9EU2Nov3Prltfw8udnUmDG/w+HI2xYWYvPG2bs5HQys3Vi+sOBpX/cTGOpj4xRyUz8SmF/i6PRHDZ0qyCUUnP7QI5+Y8zEFJqbQhQMS6JgeBKVpT7SMu3UVgUo3tnERUe9x7W3jcWVZOWlJ3eyfVM9AMmpNp5bdCqTj8pk0xoPKz+rYtKMTGaaWcm6YuPqWl55ehdpGXa+cv1ICoYl4/eFWfV5NdPnZOFO6vxj2bGlnq3r67DaLJx89mCcroRZkpKQNJb7WPagkfNh0hVDsTr189Jo4qVbBWHmf7gCmAq0jM2VUupHvSlYX5Gb7+KdDfP4YnElgwpag7Y1NYb42x0b2Lm1gb/e3hprMD3TjsNppbLMxzfO+IjkFBtlxd7o+dETU/nR3dP44qMKdm1r5NyvDuX8K4fhdFr5+J1SHv/jZr74qCJa/+G7NnLtbWP54uNKNqysJSXNxukXDmHK7GTOuGgQFcUenvnbNpJTbaxfUcuyxZWtsmQ5mHdZIcefMYihI5OZfFQm1g4yoHWFUuqIddhGwhE+/vU6gs1hcqemUzAnu79F0mgOK8SYSNRFBZG/A9+hXdpRpVTCvYrNmjVLLV++PO76LWE5bJYU3nxxTxsFAebMl2U1fPFxJRaLUDAsiVPOGYzDaeWhOzewaY0RpiE13c6w0cns3t5EY11wv/vMODab868Yxu++vwoAp8vC9GOyaW4KsWFFLS0fgdNlwe+LRNulZ9nxNoUJ+FvLXElWho1KpqEuSOleb5v7zD4ph3+9Nxeny4qnxs/W9XU0N4Zobgpjd1g4ed5gXn16F/+8fwsXfn045SVeXv/3bm78vwn8v19OIhxW2GJiE61bXkNttZ+Tz25NptPyzBJtoVx7uTw7G/nvhZ9Qub4OBE68fTL5R/d9vMeWUB6JtrgqUeWCxJUtkeVqrvBzzLVTDnShXKdviPEoiFJgIXAlxgK5i4ElSqk7DkSS3uRQK4iuCPjDfL6oguw8F+OmpmO3WwgFI7zzv2Lee6WYYaNTmDQzg0VvlVJf26o0Tjp7MOd+rZCsHGMwtn1zPf95ZDspqTa+fstYQsEwq5dWs/yTCkp2+QCYfkwWg4a4sdktnHz2IDLNtnuKGvlySSV7ihrZu6OR5qYwM47NJhiIsHFVq+JpISvXSU2lv8P/Z3Chm7JiL7mDXYwYlwoKli0xRis3/2oS379zCiKSUApi7TM7qd7SQO6kNAafmYLFZonK9d+LlrBt/j7c2Q7GX1bImPPy+2WklMidCiSeXJC4siWyXP2pIHwYeaj/DnwNyAV+rZQqOBBJepO+VBBdEWu2Kd/n5U8/XUNjfYjjz8jjG7eMxWLpvqNq9nrYtLqRUMDB7BNzsVi7blO8s4k//Wwtfq+xhtFqFQYPdeNOsuFwWqgs9VFZ5kMscNJZg9m9vRGbTZhydBYL/ruXYCCy3zUdTgvBQASl4IIrh3HsaXk89/etHHNqNj/706wem7MOJU0VPv466PXo8cRrCph7/2SysrIo+byKfx23EKvTwul/mU76sP5bFJfInQoknlyQuLIlsly9pSDimcVUZtYrw8gs5wDqD0SKgULsm+qgAje/+utMtm2oY8ax2XEph5ZrTJqZijvOL2PhyGR++LspLFtSydCRyUyakUFaZussq0hYsfzTShwOw7wVK+NJZw9m3+4mho5OocETpHRvM431QcZMTqd4ZyNP3beVN57fwxtmmOyNq+op2uDlR3dPZfJR/bPorG53EwCOVBvB5jCbntnH6AsHkXlpJh/+zJiRPfyU3H5VDhrN4U48CuJXQBXwI+CvGHGZftCLMh1xZOY4mXNK7yemGTE2lRFjUzs8Z7EKc07uWIbUdDvjp2UAkJRsY9CQ1pHU4CFuho9O4aWndlK2t5mJRyXzxYe1LHm3jCXvljF+ajpzTsll1dJqRk9M42f3Tid38KHP0BZoCuHZ2Uju5HREhHpzYkD6iGQGTc9g/b9388F31lOz3Meejyuxp9iYcPnQbq6q0Wi6Ip5prv+OOYwr1ajmyCI3383/++UkALy+Ok6al8Oi+R5WfFrFlnV1bFlnDL3Xr6jlwzf2ceHXh3PuV4cy68ScgzJDLf7telY8XMTYCwooequUpnIfVy86leFz82goMRSEK8PO+MsKKVlRRe2mJj67x4jWOvmKYXpBnEZzkMQzzfXDDoqVUur0XpBHcxiQnefgG7eM5YqbRrNqaRU7NjcwdFQyX3xUyZa1dTz39yKe+3sRuYNdzPtKIV/79mgmmCOUeAkHwnx5/1b8dUHWPLUzWr7jvTJDQRQ3A+DKcGCxCkf/fDjL796FZ0sz+bMzGXNefmeX1mg0cRKPiWluB2VH7MpqTfzY7RbmnJwXNV0df/og9hQ18tnCctZ8WUNlmY9nHyrixX/s4JFXT+SUc/LZta2B/z65g7JiL35fGJtNqCj1kZXj5M/PHhNdJLhzYQX+uiBJeU7ypqVTW9RI3a5mKjcYo5WWEYQ72wGAzWVl9q9H0rxVkTMxFenGqa/RaLonHgWRG7OfCdxBTNA+jaYFEWH42FSGj03lipsURZvqef/VEtZ8UcNNFy5hwvQMNq32EA53/H6RPyyJX90/E4DNZlrQ/NlZHHXTaCrX1/HRL9ZRs7UBpRS1uwwntTun1RFvdVgomJVYM0w0msOZeBRE7K+5HtiCkUToJ70ikeaIQEQYOymdMRPTePHxHSx6q5T1K2oRgRnHZjF6QhquJCuhkEJFFC89uZNnHtzKMafkYrdZWPXsbizAkOOM1c9pZjwsz64mrjr5Q4Z+UkU6kJSn42FpNL1FPAqiiv1NSlt6QRbNEYiIcMVNozn9wgIqy32kZzooGJa039TYqnIfH75Ryncv+ZQxwEmAB/jLf3bzy+HJpKfbsbithL1hNnxSxXizXWljiEd/v4mvfC2LIUO0U1qjOZT0NJprGNgF/Lm3BNIcmeTmu8nN73z66yXXjMBiFbZ9Vs2xlX5QsMUmbNzYwA9+vJpRo1JI94YZDEzNd2Ev9REE/vDgNjyeIDt3N3DPH8eSqnWERnPIGPDRXDWJgcNp5fJvjmJLmY/aCj9JI5K46IJ8al4spqzcT1m5n+OBwcAJo5OpLfXRBHg8RhiT0tIA99+3m8svFQLBCMlJViZOSMOqndUazQETzzTXrtJ+KqXUDYdQHs0AJtQUwrO8FgRyzxqEK93Bd64dzrsfVmK3CpMt4Pu8hsZNDYCRxNxigQvOGMRbC8tZvbqR1atjIu+m2bjp26OZe0puJ3fUaDRdEY+J6TraRXJtt9+pghCRecADgBV4Qil1T7vz9wOnmodJQJ5SKsM8FwbWmef2KKUujENWzWGMZ2UtKqRw5btw5Rm2IpfTykXnDAagaXsjxZ/XEKwOAJAz2Mnls7OYOT2DwpHCl8sbKN4TwuWy4KkLUlsf4k9/2UJllZ/zzhnMl8tq8XgCzJs3GJfOC6HRdEs8CuLPwBzgTsCCEXpjGd2sqhYRK/AwcCZQDCwTkflKqY0tdZRSP4ipfyswM+YSXqXUjPj+Dc2RQO0XRvBE94iO4yc5B7sQm6BChkssN9/N4OkZAOTl2jn/nCycpABGwMTFS2t4+8MK/vn0Lv71zK5odNtXXivhxz8cx7SpGb36/2g0hzvxKIhrgN8qpT4EEJFxwP8ppbqb5joHKFJK7TDbvQBcBGzspP6VwG/ikroTQqFQNEJrPHg8HgBslgCBUANeX+Bgbn9I8QUa+1uEDjkQuZRShGpD2LM6jzSpQora5dUAOCda8dPBfZIh97JsKl6sMo7TItF6AZrb1hU49ngnqRnZfPxJPeUVQbIyrFhtQmVVgDvu3Mhdd48mI8NGUpIVm03w+SI4HILFIjQ1hbHbBYfj4CLWNiXo55iockHiypbIcnnDfmpqa3AEeh7NtavQ/fEoCC/wBxE5FsO0dCFQHUe7IcDemONi4JiOKorIcGAkEBvWwyUiy4EQcI9S6rVO2t4I3AhQWKjzDSciNR/WUvxYCUO+VUDOWR1ndWvc0EikOYIty4Yr39HptZLGuhn09VzqVzSSMqX7oICTJyUxeVISTU1hXC4LIvDyq9Vs2OTl/366jUgEnE4L2dl29u3zM2NGCl+/Op/bf72dggInd941GptNO7o1A5N4FMS3gH8D3zCPy8yyQ8kVwMtKqXBM2XClVImIjAI+FJF1Sqnt7RsqpR4HHgcjH8SBJLKxWVJw2Bpxuw59FNKDJd5w331NT+Sq/dCIpVTxvyoKzxvZYRiM4sXG4vzkkak4peOItC04R6aQMTKn43OmiWm/8hir1RUXJvN43W727vNhtwl+f4R9+4xESqtXN7J16w683gjbt3t57616vva1YW2uFQ4r1q6rY9rU9LhnSSVaDoEWElUuSFzZElEusdrJysw60HwQnRLPNNeF5hv+BLNos1IqHltMCRAbb7nQLOuIK4Cb2923xNzuEJGPMPwT+ykITWLjK/PRtM0YmgdrA9R8UU328W0792BdkJql1SCQPrP3f3x2u4Wbrh1BQ32QtDQ7jU0hqmuCNDWHeO6VEpqbwyS5rTR7wzz7/B5efq2EyZPS+OmPxpOUZOPJf+7ktfn7OPP0PH7wvXG9Lq9G0190aWAVc7mrqRDyMRzOp8R57WXAWBEZKSIODCUwv4N7TMCI8bQ0pixTRJzmfg5wAp37LjQJTM1nhr9A7Mabdulr+78jVH5Yjgop3EPduPooRLfNKmRmOrBahfQ0O6NGJDF1UhqnnZiD22Xh8vPyOXpaOpEINDWF+XJZLb/6zQZWr/Ew/819ALy/sILPPjOsrV98Wc37C8uJRBT33LuZa29YxnsfGMftaWoKdRiPatfuJqqrO04Jq9H0B52OIERkIcY01jNE5AZMM4557jdKqd91dWGlVEhEbgHexZjm+pRSaoOI3AksV0q1KIsrgBdU29ynE4HHRCSCocTuiZ39pDl8qP7EUBDZJ+VSvbiSxs0NePc14y5IitapWlQBQKqZDKg/OWtuLmfNNdZNjBuXwvGzMgmFFf/+XwmbtzTwi18b6ywy0+3U1gX5031bmPpuGitWeQBYuLCCteuNiLN/fXAbY992880bhjB9kjEyWvRxBfc/sI1pU9O5647J7C32snxFLZ8trWbjpnqyshz845GjcbvbTsMNhSKsWVvHuLEppKYeWjOCRtMZXZmYpgAt6xa+Y27vAk4Gvg10qSAAlFILgAXtym5vd3xHB+0+A6Z2d31NYhP0BGgqakRsQvqMdAJVfurX1lH2Zikjbxxt1KkP0ryrGbEKqZO79j30NVaLMMTMU/7da4fz+ttlbNnRRJLbyreuGsb7iytZvaGeFas8WK1i+CZM5XD01DQ2b29i2zYvP/9ZEdOmVpKaauOzpdUoBStXeXj8iZ288dY+IjHpwGtqArz0v2KuuXp4tKy5OcTd92xm1WoPLqeFM04fxBmn5zF2TEq/K1TNkU1XCiIdqBaRdAz7/x6l1B0ici3waJ9IpzmsadhsrHh25jmxJdlIn5FB/do6qj6uZMQNoxCr0LDRSG/uHOTE5opnzkT/kJXp4PqrhlHjCaLCiuxsB1dcMoTjZmWxfHUt0yamUVrpZ8HCCqaOT+WyC4YQCIZ5e9E+VqxuZK2ZdQ9gaL6LvaU+Xn/DMFWNGZ7EuFHJZGY6eO6VEl59rYTRo5KZMiUdt8vKL369nq3bGrHbBJ8/wpsLSnlzQSmnnJzDD24bd9BTcTWazujqF7kLIw/15RhmnnfM8mHEN81VM8Bp3GIqCNOv4B7qxp7lIFgToGpxJbmn5lFvJgBydhHIL5HIymhr3hk+1M3woYbsY8ekMH1yOilJFiwWYxX4ufMyOXVuOps3hImEFAWDnRQUJHH/o9uprQsyapiba742NNrJj1/jYcv2Ju6+ZzPJyVYmTUxj67ZG0tNsfOPSQsQmfLG8llUb6vh4cRV793o5+aQcTj81j+zstqHPN2+uJz/fTXr6kWuSUkrx9LO7yc5ycMH5Bf0tzhFHVwri18CzGKaeKuAvZvkVwOe9LJfmCKBhizE6cJlmGhEh4+hMKt8vZ/sDWwnVB2nYYNRxFx4eCqI70tP2/0m5XRaOPTqtTdlVlw1h5Zo6Tjkuu80I4GsXD+HDJZVs39VMaYWfZctrsVrg0nn5FJrP6NLz8zl2dib//M9eduxsYsfOJl76XzGnnpLHylW1FBa6SU62seijSnJzHPztrzPZW+wlI8MeNZl1RHW1H7EIWZn7r0OJRBTLV9QybGgSgwcbCr+mJsCu3U3MnJHRb6auNWvr+O/LxYhghFwZmtR9I03cdKoglFIvmfmoRwGblFKNImIDrsJYC6HRdIoKq+j01qQRrZ1S5uxMgnUBPF/WsvtJM9e0BZJHdRxe40hlaIGboR101kluK+efNZiIUryzsIKly2s59bhsxo9ru76jYJCLH353FBs2N7BijYede728ucBYS7Kv1BetV1kV4Du3rMTjCRrJmqZnsK2okZQUC+eel8MZp7hZsbKWl18pYfeeZkRg+tR0LrqwgFlHZ1FfHyQlxcbDj2znvQ/KEYETjsvmuGOzefzJndTVBfnmtSO48IICSvZ5yR/swuVqdbCXlvlYt87DSSfmRh3v69bXEQhEOPqozLifVzisOlxzsuAdoytSCp5+djejx6RQWeVn8sQ0TjoxB7tdm98OBmk7eejwZtasWWr58uVx128Jy2GzpPDmi3sY1MXbVV/j9Rmml0RbKBevXE3bG1n3g9XY0m2MumXMfm+Ydes8lM0vBQWOPCcjbxx1UHK1hNzobKFcf3GwckUiChGhqxd0pRSfL69lx+5mpoxLZfc+L6XlPmZPz+SND8rw+iLY7WJk7+vi5+6wC+EI0Sm4FgtEIq1bm1WIKNXGqd5SLyvLQVVVABFjZbrNKuTludi1u4lIBMaNSeF3d06hoSHITTevJBRS3P3bKcycmQGAzx9m584mysp9NDQ1MX1GKsMLcqirC3LfA1vZsLGeW28ewykntUbmrakNcO03l0WnErf/32ZMT+eOX08+ZD6aBvO7n2gL5Rp8dTRX+Dnm2ikHulCu029X4noFNYctEX+Yig/KAcP/0JH5IX1qBigof7uMlPGJNXspkbBYujfdiAjHzc7iuNlGFIHp01o7sNxcB8tXezhmZgYOl5Ut2xoozHdTWd/A8pWN7C0J4nZaOPmYLObMyjKUzYpaPv2ihsbmME6HBX/AUDCXn5vPsGFJLF5azZoN9Ywc6iYl1cbnKzxUVQVITrLi9YXx+QwN0rizCYsYJratRY389OdrycpyEDKDLf7pL5v5630zWL68ln88tZNAoFXzJCdbOf+8Jt7/oJyaGmNd7h/v3cKCt8vIzXUwY1oGSz6tIhxWjB2RRHKKjdXr68lMtzN1fCrL13pYvaaOu+7eyO2/mtTlSMLjCVCyz8vIEckkJXXcJQaDEZZ+5iE7x87R0xJLQfQmegSBHkH0hO7kCjUEWfejNfjLDDNH3lmDyJzTefgTpdQhsV8fqSOI3qJFLms4GYH9zDcRpQj4FU6nBb8/TDikSE7Zv/MMhRUfLqnEbbdyzOxMrFYhEIgQDCkqq/ykJNlwuiw8/uxuaszkTnabkJvtYF+5H5HWN/+cLAc5mXYamgOUlAaj9xgy2MnYkSl8/Hn1fqMEq1W4+pIhjB6VzNqN9YwdlUx6mp2yCh+PPbMbry/CnFmZfPc7o6mrC7JlawNTp6QzbGgSq1Z7eP7FvWzcZPjBsrIcXHJRAe+8V05dXZC8XCffu3UstbUB/vb3IqqrA1it8NvbJ2OxCB5PkCFD3Iwd07+fbW+OILSCQCuIntCdXEX3b6VqUQX2DDuZx2SRcXQmEsdb8MGS6B3xQJer2Rvm+VdK2LaziZPnZHHicVm8tqCMzdsbsYhw9ik5nHBsDhYLeFUDSz6pZ/euEDMnpTN9Whp2u4Wa2gD7ynxUVPop2tVMaoqVU47JpmBIx7/bfWU+Hn92Nz5/W5uYiDFtudocmdjtgttppb4xtN81kpKs+HxhIhFISrLQ3BzZr85llwzhhutHduonaUEpxZq1dQwpcJOb66SuLkhyspVgUPHQ34sYPjyJr35lKCtX1ZKWZmfM6Pg+m4raGlZ/7uHGPx3b9wpCRMYDPwZGYKyIBiOT3OkHIklvohVE79OVXJ6VtWy+YwNiEwqvGkrSsL5zPOuOuGf0h1xKKWpqg2RmOLCYFh+vL0wgECE9rbVjO5SylZT5mP9OGdW1AcMvku2gaHczkQikpdiYOTmNE4/Nxp1k5a33ytm6o5FjpmcwbWoarywoY+v2JgDmTE/n7HkpvPpGDes3eklPtTEox8G2Xc0oBcfMyWL5iloG5TkZPy4Vu91CszdMbW2A6mo/Rx+dRSSsePvdMhwOC1Mnp7FqjYfCIUnk57v44kujLzr+2Gw++7wam024955pjB/X1vwaDit272mioMBNQ32I517Yw8eLKwkEIry/YR4jJh5Qf3FQCmITML5dsVJKJVxKLq0gep+u5Cr6yxaqPq4kY1Ymg+YN7lO5dEfcMxJVLuh92erqg1RVBxk+1IXN1rlvIhiK8M7CClLcVk45MZegpZGIUngqDFOY3WFh0adVvLuoMu57tzj822O1QLhdeWaGnTvvmEx6up1Vqz3sK/Xx8ccVlJX7SUm2Eo6A12sEwC4cZOPe509k9qmD4pYlhoNyUmcB9wN/wsjNoNF0iL/SCDSXNFzPRdckLulp9jYjls6w2yxccHbbFx2LCPmDWhcknnJ8NrWeILWeIKcenw0CZeV+QuEIDruVlGQrTqeFDxZXUVUT4KKzB+FOsrFteyOTx6Xy3uJKdu5t5oIzBlFa5efLVR5Omp3F7n1edhd7ufX7q/eTy+W00NhkKIaxI5OYe3IKySrMmIlp+9U9WOJREM8A44AUoPaQS6A5rIn4w/gr/bgLk/BXGQrClnHkrtzVaGKxiHDpefltykZ1kDJ37OiUNhMyxprrfr519TB8vggulxUROGtuHkluYzbYe4sqWbHWQ0TByKFucrOdDBnkYvoUI6xLU2OI8WNTCVoa8ext3u+eh4J4FMSPMKK6nh9TpuJsqznC2f2vXZS/Vcqku6cSMBWEI6vzjHAazUClo9l6ItImcm9ykrGf5LZy8bmDOf+sPMIhhdPV1qI/rBPH/KEmnk5+MYZC0Gj2o26NB4DKRRUQAYvbitWZcO4pjeawxGazYOvHV/F4MsrN7QM5NIchEX8Y3z4vAPVmmGtbB3PlNRrN4Um3v2Yzq9wVGEH7WtJ9KaXUj3pTME3i07zXC+bMi5aFcbYUPXrQaI4U4nndexgjYZCidTqUwvBNaAYwzTsb9yuzJusRhEZzpBBPFKtLgP+Y+98DFmFkltMMcJp37T9zwq7TYWo0RwzxKIhMYIm5Xwq8DNzYaxJpDhuadzXtV2ZL1yMIjeZIIR4FUYZhiioDnsBIHBRX/FwRmSciW0SkSER+1sH560SkUkRWm3/fijl3rYhsM/+uje/f0fQF3l1eSv67l6adhoJwxST7sR/B2cs0moFGPK97v8LIKPcj4K+AF/hBd41ExIrhvzgTKAaWich8pdTGdlVfVErd0q5tFvAbYBaGv2OF2VYv1OtnIsEIO/+0m2CVEW3T4rSQPCoZX7Exm0mvgdBojhzimeb6bwARyQCGK6X8cV57DlCklNphtn8BuAhoryA64mzgfaVUjdn2fWAe8HxXjUKhUDS+Ujx4PB4AbJYAgVADXl8g7ra9jS+wvwM4ESh/r4xgVRBxCiqgcA53YMkxl8lYIJLux0+w64v0AgF6ZyXpwaLl6jmJKlsiyxW2efHU1aIc3h63z8rqPBx/t6YiERkhIsswRhEnicjHInJnHPcdAuyNOS42y9pzmYisFZGXRWRoD9siIjeKyHIRWV5dXR2HWJoDJRKIUP26B4CM09IZ9tMh5F2ajT3HeM+wJluQLoKfaTSaw4t4TEyPYnTOgjHrfTHGuojbD8H93wCeV0r5ReQm4GngtJ5cQCn1OPA4GNFcu9KGnWGzpOCwNeJ2JU401xYSKZpr5aflhD1h7Lk2cqYPwmIqA2euImduBGuytd+jg/b3/TtDy9VzElW2RJTLGoKM9Ewysw5tHxbP697xwEMxx9uBwjjalQBDY44LzbIoSqnqGJPVE8DR8bbV9D1VHxlhjZOmJEWVAxjxZLJPzCFjZvxJ6DUaTeITj4KoAqaY+3kYo4d9cbRbBowVkZEi4jDbzY+tICKxYRAvBDaZ++8CZ4lIpohkAmeZZZp+IlDtp26tB6yQOlWH89ZoBgLxmJj+Adxt7j9nbvebstoepVRIRG7B6NitwFNKqQ2m/2K5Umo+cJuIXIiRZ6IGuM5sWyMid2EoGYA7WxzWmv6hanElKHANd2KLI5a+RqM5/IlnFtMfRGQfcJ5Z9KZS6pl4Lq6UWgAsaFd2e8z+z4Gfd9L2KeCpeO6j6V18+7yUvmJY+JIm6NGDRjNQiGvZq1LqaQwHsmaAEWoMsen29QTrgriGuEmZmniOfI1G0zt0qiBEJNxFO6WU0jEVBgA1n1Xhr/DjyHFQ8JUhROzxLoPRaDSHO1118oKxinkf4OkTaTQJR+PWBgCSx6ViT7XjRysIjWag0NUspn8CTUAOsA74oVJqastfn0in6ReUUnhW1BJuDtG4zVjR7S5wddNKo9EcaXSqIJRSNwD5wP/DWJPwjojsEpF5fSWcpn+oeLeczb/dQNFft9G8uwkE3EO1c1qjGWh0uQ5CKdUE7AB2AgGM0URqH8ilOUT4ynyEmkJx11dKUf6Wscyl9vNqiIA904FNJwLSaAYcnSoIEfmliGwDPgTGALcC+Uqpl/pKOM3BEaj2s+bmFay5eSW+0viCeDVubqB5d9ugZM48Z2+Ip9FoEpyuXgvvwnBS78BYTX0hcKGRohqllLqo98XT9JRAtZ9Nd2wgc04WScOSUUFFsCbAhp+vY/I903AN7tiX0LitgaL7thL2GpPX7Jl2grVGVFbnIK0gNJqBSHd2AwFGm3+xqN4RR3OwFD+/B+/uZvylPnJOzQNArEKwJsDGn69l8h+n4cxzUTq/BIvdwqBzjGgnu/+5E1+JOcqwwOAL8il5sZiIP4K7UPsfNJqBSFcKYmSfSaHpkMq3qrCG6xny1aHdVwa8xc1UfFAOGKG5qxZVAJB3zmDqVtTiK/Wx5XcbGX/7ZHY/sRMEMmZn4i/307C+HovTQt5Zg7AmWXEPTaLgskK8xc0kDdcKQqMZiHSqIJRSu/tSEE1bIqEI+54thYjRwdtTu49/VPLiXoiA2AQVUkQCERBIGZdC6oRUdjxYRPOuZqo+MhQHCqo/rqRubR0AqZPTSJ+eEb1e8qhkkkcl98a/p9FoDgN0dpc4UEqx55ldlPyv+IDaV39axcZfr2PL7zdSuyy+mIPByqCRfQPw7u0+k5WKKDwrjYys2SfnRMsdOU5sSTasLitJI4yRQOmrrZHT9722j7pVHsQuZB2XHe+/pNFoBgBaQcRB884m9r1czN5ndhlv5T0g7A+z8+9F1K+po/bzGrY/sBUV7t6F4y9vTX/qLe5+BpJ3bzOhhhDWZCsZszMRmwDgjHFKJ40wRgOhhtZpr6E6wxGdPiMDR6bOJ63RaFrRCiIOqpdUGTvKsPP3hMqFFYQaQjhyHNjS7YTqQ9E3/a4IxCgI377uFUT9esNM5CpwY7Vbo6Yh95AYBTGy1VwkVsFt+hasKTayT8xBo9FoYtEKohuUUlQvqYweN+9uAgwfQfVnVUT8ncc0VGFF6WuGOSf96EzSpxvpQyveK2tTL+IPEwm2HZnEKgh/ma9bOes31AOGggAYdG4+g87LJ216a8pSR44Da7LV2M9zkn1CDtZkKzmn5OqFcBqNZj+0guiGxq2N+CtaA9R5S4zOuvKDCrbds5kdj2zfr42v3IcKK2qWVuMv82FLt5M+PZ20KUZn7VlZG11vEPaGWf3dFWy6fX2ba8SamPzlnSsIpRTh5hAN5ggiaaQxKrCl2MiYmYHVbo3WFZHoKMKV7yJ5VDJjfjCOjJkZcT8PjUYzcNCvjd1Q85lhXrI4LEQCEXwlhompqciIclr9cSXDrh0Rtd83bKxnw8/WkjE7k6DHsO+nTU3H6rBizbLiKnTjK/ZStaSSQWcNpn5jHYGqgPFX68eRaSxKazOCqGrdb8+OB7dRudCYlWR1W3Hldx1UL+fkXMQqZB2jHdIajaZr9AiiGxpM003qJCMElc8093j3Gn6BWDMSQLWpUDzLamna1ojFbSVzdmb0fNpUYxTR0qnXr6uLnqtfV0fN59XULq9poyBC9UEioVYTVPUnVRQ/v4ewP0xVi38ESB6bgsXa9UfqyHKQf0EBjmztkNZoNF2jRxBdEPaHadrRCAJp0zKoW11HoNIwN8U6qyveLaPwimFY3VbqVnnaXCNtUlob+37axFQq3i2jcXM9gdpA1LkMUL24qnUarAJxClanlVB9CF+Zj6TCJHzlPoru24IKKULeECoQwZ7tYPgNI7DYtb7XaDSHjl7tUURknohsEZEiEflZB+d/KCIbRWStiCwUkeEx58Iistr8m9+bcnZG07ZGVEhhz3LgHuIGgVB9CH+ln1BDCLELrgIX4eYwJf/dg7/Kj3dvM2IXMuZk4shzknlsVptrWpNspIxJAQVl80toKmqMnqv9ssYIYmLOgrWlW7GbpquWtRDF/96NChkVyuYbUVfdhW6sDitmnCyNRqM5JPTaCEJErMDDwJlAMbBMROYrpTbGVFsFzFJKNYvId4E/AV8zz3mVUjN6S754aNhozgwa7EKsgj3DCGBX+0U1YITBzj1jEHuf2U3p/H1YnIZD2FXgZtBZgzu9btrUdBq3NlL6+j6IYEx/NdcjxGJNs+FItuPdDb4SL/Xr66j6uNJQ6xGiC+la1jdoNBrNoaQ3TUxzgCKl1A4AEXkBuAiIKgil1KKY+p8DVx/MDUOhEDU18a1UBvB4PADYLAECoQa8vrbOYM8G41q2QsFPI7YsK8HaIFWfG/GOrJmCdViEpAlumjd7Kf7PHgDsw6z4aaQz7OME51AH/r3G/VxjHPh2KEK1ISxuC84hDrxFPiwZESTFaFP29j5KXt4LQMr0JIKVIfzFAbCCfQxd3u9QEqBn60D6Ci1Xz0hUuSBxZUtkucI2L566WpQjvrD+sWRlZXV6rjdNTEOAvTHHxWZZZ9wAvB1z7BKR5SLyuYhc3FkjEbnRrLe8urr6oASORUUUzVuNNQ/u4cbMIlumoU+bNhpfFFuWcZx1VgaOAjNWkkDSWHeX1xarMPjqPNKOScGaaiV5khtnoWFKco9xkX1eFsmznaTMdpE6Ixlbpo1gZZBIcwTXKCdZZ2WQPNUYNTgLHFjd1q5up9FoNAdEQjipReRqYBZwSkzxcKVUiYiMAj4UkXVKqf0WHSilHgceB5g1a5bqSht2hs2SgsPWiNvV2rH7K32EmyJY3VZScrMQhPSx0LC8MeoDSMpKwUkKzjQY+c0MvPu8hJvCJOcnI3TjD7BC/pmphgEOSErxU62qyDw+G1eqC9vZRqfvJIWRN6ZS/UkVoYaQEW3VbsU1Mw1r0I57qBsnfR9t1UlKn98zHrRcPSNR5YLElS0R5bKGICM9k8ysrl9Oe0pvKogSjFzWLRSaZW0QkTOAXwKnKKWiK9KUUiXmdoeIfATMBPZfldZLBKoN84811RZ1/iaPTiF1clp06qszr+2aA3fBgX84jmwn+Rd3PMCy2C3kmrkdWhCLkHWsXsug0Wh6j940MS0DxorISBFxAFcAbWYjichM4DHgQqVURUx5pog4zf0c4ARifBd9QYuCsCW3Nd8MmjcYW6oNa5JVp+LUaDRHNL02glBKhUTkFuBdwAo8pZTaICJ3AsuVUvOBe4EU4CXzLX2PUupCYCLwmIhEMJTYPe1mP/U6gWpjMGNtF6PI6rYy8jujUGGl1x1oNJojml71QSilFgAL2pXdHrN/RiftPgOm9qZs3RE1MXUQxK5lOqtGo9EcyehX4E4I1BgKwp6aEH58jUaj6XO0guiEFhOTLa37VJ8ajUZzJKIVRCcETROTPV0rCI1GMzDRCqIDlFJRH4Q9QysIjUYzMNEKogPCTWEigQhiF6xJ2iGt0WgGJlpBdEDsFFcdIVWj0QxUtILogNZFcnoGk0ajGbgMeAVRv8fLZ7/dQN1CI0JrxXtlVH1kLOpuv4pao9FoBhID/hW5YY+XFX/dijXDTsPRWex4qCh6rqNFchqNRjNQGPAjiILjM0nOdxH2BNn5j7axALWC0Gg0A5kBryDEIoy+sACA5u1G/of0ozNwFbhInZB4YX01Go2mr9CvyMDYi4aw9rEdANiz7OSdNQiLdcDrTo1GM8DRvSCQNTEVa7qhK5PHpGjloNFoNOgRBAAiQvLRWUT2eck8pucZ6TQajeZIRCsIE/tgFznH6wxtGo1G04K2pWg0Go2mQ7SC0Gg0Gk2HaAWh0Wg0mg7RCkKj0Wg0HdKrCkJE5onIFhEpEpGfdXDeKSIvmue/EJERMed+bpZvEZGze1NOjUaj0exPrykIEbECDwPnAJOAK0VkUrtqNwC1SqkxwP3AH822k4ArgMnAPODv5vU0Go1G00f05jTXOUCRUmoHgIi8AFwEbIypcxFwh7n/MvCQGAkYLgJeUEr5gZ0iUmReb2lXNwyFQtTU1MQtoMfjASDSZCdobaY5EI67bW8TlGYAwirSz5K0RcvVM7RcPSdRZUtkuUIWP566WpTD2+P2WVmdr/3qTQUxBNgbc1wMHNNZHaVUSETqgGyz/PN2bYd0dBMRuRG4EaCwsPCABLW7bKRkO/DWBQ+ofW8QsoaMbThxZAItV0/RcvWcRJUtkeVyZzuwOw+9QeiwXyinlHoceBxg1qxZqitt2BlZWVlcf3/+oRbtoGgZCR3I/9ObaLl6hpar5ySqbANRrt50UpcAQ2OOC82yDuuIiA1IB6rjbKvRaDSaXqQ3FcQyYKyIjBQRB4bTeX67OvOBa839rwAfKqWUWX6FOctpJDAW+LIXZdVoNBpNO3rNxGT6FG4B3gWswFNKqQ0iciewXCk1H3gSeNZ0QtdgKBHMev/FcGiHgJuVUonjQdZoNJoBQK/6IJRSC4AF7cpuj9n3AZd30vZu4O7elE+j0Wg0naNXUms0Go2mQ7SC0Gg0Gk2HaAWh0Wg0mg7RCkKj0Wg0HSLGrNIjAxGpBHb3sFkOUNUL4hwKElU2LVfP0HL1nESV7UiUq0opNa+jE0eUgjgQRGS5UmpWf8vREYkqm5arZ2i5ek6iyjbQ5NImJo1Go9F0iFYQGo1Go+kQrSDMQH8JSqLKpuXqGVqunpOosg0ouQa8D0Kj0Wg0HaNHEBqNRqPpEK0gNBqNRtMhA1pBiMg8EdkiIkUi8rN+lGOoiCwSkY0iskFEvmeW3yEiJSKy2vw7tx9k2yUi68z7LzfLskTkfRHZZm4z+1im8THPZLWI1IvI9/vreYnIUyJSISLrY8o6fEZi8KD5nVsrIkf1sVz3ishm896vikiGWT5CRLwxz+7RPpar089ORH5uPq8tInJ2H8v1YoxMu0RktVnel8+rs/6h979jSqkB+YcRgnw7MApwAGuASf0kSz5wlLmfCmwFJmHk6/5xPz+nXUBOu7I/AT8z938G/LGfP8cyYHh/PS/gZOAoYH13zwg4F3gbEOBY4Is+lusswGbu/zFGrhGx9frheXX42Zm/gzWAExhp/matfSVXu/N/AW7vh+fVWf/Q69+xgTyCmAMUKaV2KKUCwAvARf0hiFKqVCm10txvADbRSQ7uBOEi4Glz/2ng4v4ThdOB7Uqpnq6gP2QopRZj5DOJpbNndBHwjDL4HMgQkV7Jd9uRXEqp95RSIfPwc4xsjX1KJ8+rMy4CXlBK+ZVSO4EijN9un8olIgJ8FXi+N+7dFV30D73+HRvICmIIsDfmuJgE6JRFZAQwE/jCLLrFHCY+1demHBMFvCciK0TkRrNskFKq1NwvAwb1g1wtXEHbH21/P68WOntGifS9+ybGm2YLI0VklYh8LCIn9YM8HX12ifK8TgLKlVLbYsr6/Hm16x96/Ts2kBVEwiEiKcD/gO8rpeqBR4DRwAygFGOI29ecqJQ6CjgHuFlETo49qYwxbb/MlRYjle2FwEtmUSI8r/3oz2fUGSLyS4xsjc+ZRaXAMKXUTOCHwH9EJK0PRUrIzy6GK2n7ItLnz6uD/iFKb33HBrKCKAGGxhwXmmX9gojYMT7855RSrwAopcqVUmGlVAT4B700tO4KpVSJua0AXjVlKG8Zsprbir6Wy+QcYKVSqtyUsd+fVwydPaN+/96JyHXA+cDXzY4F04RTbe6vwLD1j+srmbr47BLhedmAS4EXW8r6+nl11D/QB9+xgawglgFjRWSk+SZ6BTC/PwQx7ZtPApuUUvfFlMfaDS8B1rdv28tyJYtIass+hoNzPcZzutasdi3wel/KFUObt7r+fl7t6OwZzQeuMWeaHAvUxZgJeh0RmQf8FLhQKdUcU54rIlZzfxQwFtjRh3J19tnNB64QEaeIjDTl+rKv5DI5A9islCpuKejL59VZ/0BffMf6wgufqH8Y3v6tGNr/l/0ox4kYw8O1wGrz71zgWWCdWT4fyO9juUZhzCBZA2xoeUZANrAQ2AZ8AGT1wzNLBqqB9JiyfnleGEqqFAhi2Htv6OwZYcwsedj8zq0DZvWxXEUY9umW79mjZt3LzM94NbASuKCP5er0swN+aT6vLcA5fSmXWf4v4Dvt6vbl8+qsf+j175gOtaHRaDSaDhnIJiaNRqPRdIFWEBqNRqPpEK0gNBqNRtMhWkFoNBqNpkO0gtBoNBpNh2gFodF0ghmxU7X78/TCfe4wr/2VQ31tjeZgsPW3ABrNYcAqjMiZAIH+FESj6Uv0CEKj6Z5KjIVIHwALReQ6843/32YugCoR+XFLZRH5thmjv0lEvhSRE81yh4j8QUR2m7kEFre7z6li5GqoFJHLzTYnmAHsfGZ5n0cT1QxctILQaLrnLAwlUUnbsCKnYgSZKwPuFZHpInIaRgL5SowgbsOA+SKSjRGz/2cYK3BvwViBG8vp5vXSgXvMsp9irGi/GbgTqDrU/5xG0xnaxKTRdM8XwK/M/Vpgqrn/lFLqMREJAU8Ap2AoBIDfKKXeF5FhwC8wErdcgBEy4WvKiOvfnvuUUo+LyHcxYvuAEUbhfIzQCisxQihoNH2CHkFoNN1TpZT6wPxbEVMu7baxqHbbeGhJVhOi9bf5fxiRRLdhxCxaLmaaUI2mt9EjCI2mewpE5IqYY7u5vV5E9gC3mccfYwRQ+xHwWxEZjdGp12Jkb3sDmAW8KCIvA9OUUt/v5t4/B/wYZqm9GGk30wDPQf5PGk23aAWh0XTPTNomi/mBuV0I/D9gMPATpdQaADPz3k+B+4CNwA+UUtUicg/gBr4OnEZ8YasjwK3mPaoxciLvOej/SKOJAx3NVaPpIWbCnX9iKIU/97M4Gk2voX0QGo1Go+kQPYLQaDQaTYfoEYRGo9FoOkQrCI1Go9F0iFYQGo1Go+kQrSA0Go1G0yFaQWg0Go2mQ/4/40MN50/3NSEAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Train MSE: 0.6208\n",
      "Test MSE: 0.7135\n"
     ]
    }
   ],
   "source": [
    "from tqdm import tqdm\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import torch\n",
    "from torch import nn\n",
    "import torch.nn.functional as F\n",
    "from torch_geometric_temporal.nn.recurrent import GConvLSTM\n",
    "\n",
    "from torch_geometric_temporal.dataset import ChickenpoxDatasetLoader\n",
    "from torch_geometric_temporal.signal import temporal_signal_split\n",
    "\n",
    "loader = ChickenpoxDatasetLoader()\n",
    "\n",
    "lags = 10\n",
    "stride = 1\n",
    "epochs = 200\n",
    "\n",
    "dataset = loader.get_dataset(lags)\n",
    "\n",
    "train_dataset, test_dataset = temporal_signal_split(dataset, train_ratio=0.2)\n",
    "\n",
    "### MODEL DEFINITION\n",
    "class RecurrentGCVAE(nn.Module):\n",
    "    def __init__(self, node_features, window_length):\n",
    "        super(RecurrentGCVAE, self).__init__()\n",
    "        \n",
    "        self.n = window_length\n",
    "        hidden_state_size = 64\n",
    "        latent_dim = 16\n",
    "        \n",
    "        self.encoder_rnn = GConvLSTM(node_features, hidden_state_size, 1)\n",
    "\n",
    "        self.latent = nn.Linear(hidden_state_size, 2*latent_dim)\n",
    "\n",
    "        self.decoder_input = nn.Linear(latent_dim, hidden_state_size)\n",
    "        self.decoder_rnn = nn.LSTMCell(input_size=1, hidden_size=hidden_state_size)\n",
    "\n",
    "        self.reconstruct = nn.Linear(hidden_state_size, 1)\n",
    "\n",
    "    def forward(self, window, h, c):   \n",
    "        x = window.x\n",
    "        edge_index = window.edge_index\n",
    "        edge_weight = None\n",
    "    \n",
    "        H, C = [], []\n",
    "        \n",
    "        for t in range(self.n):\n",
    "            x_t = x[:,t].unsqueeze(0).T\n",
    "            h, c = self.encoder_rnn(x_t, edge_index, edge_weight, h, c)\n",
    "            \n",
    "            H.append(h.detach())\n",
    "            C.append(c.detach())\n",
    "            \n",
    "        mu_z, logvar_z = self.latent(h).chunk(2, dim=-1)\n",
    "\n",
    "        std_z = torch.exp(0.5 * logvar_z)\n",
    "        eps = torch.randn_like(std_z)\n",
    "\n",
    "        z = std_z * eps + mu_z\n",
    "        \n",
    "        h = self.decoder_input(z)\n",
    "        c = torch.zeros_like(h)\n",
    "        \n",
    "        x_hat = self.reconstruct(h)\n",
    "        x_rec = [x_hat]\n",
    "        for t in range(self.n-1):\n",
    "            h, c = self.decoder_rnn(x_hat, (h, c))\n",
    "            x_hat = self.reconstruct(h)\n",
    "\n",
    "            x_rec.append(x_hat)\n",
    "\n",
    "        x_rec = torch.cat(x_rec[::-1], dim=1)\n",
    "        kld = torch.mean(-0.5 * torch.sum(1 + logvar_z - (mu_z ** 2) - logvar_z.exp(), dim = 1))\n",
    "\n",
    "        return x_rec, kld, H, C\n",
    "    \n",
    "model = RecurrentGCVAE(node_features=1, window_length=lags)\n",
    "optimizer = torch.optim.Adam(model.parameters(), lr=0.01)\n",
    "kld_weight = 0.5\n",
    "\n",
    "### TRAIN\n",
    "model.train()\n",
    "\n",
    "loss_history = []\n",
    "kld_history = []\n",
    "\n",
    "for _ in tqdm(range(epochs)):\n",
    "    h, c = None, None\n",
    "    avg_rec_loss = 0\n",
    "    avg_kld = 0\n",
    "    for i, window in enumerate(train_dataset):\n",
    "        optimizer.zero_grad()\n",
    "\n",
    "        x_rec, kld, H, C = model(window, h, c)\n",
    "\n",
    "        rec_loss = torch.mean((x_rec - window.x)**2)\n",
    "        loss = kld_weight * kld + rec_loss\n",
    "\n",
    "        h = H[stride-1]\n",
    "        c = C[stride-1]\n",
    "        \n",
    "        avg_rec_loss += rec_loss.item()\n",
    "        avg_kld += kld.item()\n",
    "        \n",
    "        loss.backward()\n",
    "        optimizer.step()\n",
    "    avg_rec_loss /= i+1\n",
    "    avg_kld /= i+1\n",
    "\n",
    "    loss_history.append(avg_rec_loss)\n",
    "    kld_history.append(avg_kld)\n",
    "train_loss = avg_rec_loss\n",
    "\n",
    "### TEST \n",
    "model.eval()\n",
    "rec_loss = 0\n",
    "with torch.no_grad():\n",
    "    h, c = None, None\n",
    "    for i, window in enumerate(test_dataset):\n",
    "        x_rec, _, H, C = model(window, h, c)\n",
    "        \n",
    "        rec_loss += torch.mean((x_rec - window.x)**2).item()\n",
    "        \n",
    "    rec_loss /= i+1\n",
    "test_loss = rec_loss\n",
    "\n",
    "### RESULTS PLOT\n",
    "colors = ['#2300a8', '#8400a8'] # '#8400a8', '#00A658'\n",
    "plot_dict = {'Reconstruction': (loss_history, colors[0]), 'KL-Divergence': (kld_history, colors[1])}\n",
    "\n",
    "n = len(loss_history)\n",
    "\n",
    "# plot train and val losses and log_likelihood area under the curve\n",
    "fig, ax = plt.subplots()\n",
    "x_axis = list(range(1, n+1))\n",
    "for key, (data, color) in plot_dict.items():\n",
    "    ax.plot(x_axis, data, \n",
    "                label=key, \n",
    "                linewidth=2, \n",
    "                linestyle='-', \n",
    "                alpha=1, \n",
    "                color=color)\n",
    "    ax.fill_between(x_axis, data, \n",
    "                alpha=0.3, \n",
    "                color=color)\n",
    "\n",
    "# figure labels\n",
    "ax.set_title('Loss over time', fontweight='bold')\n",
    "ax.set_xlabel('Epochs', fontweight='bold')\n",
    "ax.set_ylabel('Mean Squared Error', fontweight='bold')\n",
    "ax.legend(loc='upper left')\n",
    "\n",
    "# remove top and right borders\n",
    "ax.spines['top'].set_visible(False)\n",
    "ax.spines['right'].set_visible(False)\n",
    "\n",
    "# adds major gridlines\n",
    "ax.grid(color='grey', linestyle='-', linewidth=0.35, alpha=0.8)\n",
    "\n",
    "# y-axis in log scale (reconstruction loss tends to start high)\n",
    "# ax.set_yscale('log')\n",
    "plt.show()\n",
    "\n",
    "print(\"Train MSE: {:.4f}\".format(train_loss))\n",
    "print(\"Test MSE: {:.4f}\".format(test_loss))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8597420e",
   "metadata": {},
   "source": [
    "### GConvLSTM encoder - Linear decoder"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "fa244ffb",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 200/200 [26:34<00:00,  7.97s/it]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEWCAYAAABrDZDcAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAABeiklEQVR4nO3dd3hUVfrA8e+Zkt4bIaEkgdB7FUUpCmJDEVFcG2J3LWvbddWfy7qWtZcVO/beRURRAQVBeg8kIYQkhIT0Xqec3x/3JkxCEhLIZAbmfJ4nT2bu3Lnz5ia57z1dSClRFEVRPJfB1QEoiqIorqUSgaIoiodTiUBRFMXDqUSgKIri4VQiUBRF8XAqESiKong4lQgU5QQjhPhNCCGFEPNcHYtyclCJQHEbQogM/QJ3katjcQdCiMn6+cho9tKXwIvA7q6PSjkZmVwdgKIoIIQwSykt7dlXSvmys+NRPIsqESgnDCHELCHERiFEhRAiUwixUAgRor/mJYR4UwhxSAhRJ4Q4IIT4Xn9NCCEe17fV6fssE0KEt/I5/kKIp4UQ+4QQlUKIbUKIq/TXEoQQdiFEkRDCrG/rrd+5F+lxmIQQfxdC7BFCVAkhdgshbnQ4/gJ9/y+FEJ8LIWqAK5rFMBlYqT9tOL7UX2tSNSSEeFd//r4Q4kchRI0Q4mc9rq/0GP4UQsQ7HH+IEOIHIUS+EKJA369XJ/yalBOQSgTKCUEIcS7wNTBM/14B3Ap8qu9yNXA9UAgsAjYDp+qvnQn8E7Dpr60ChgKBrXzcO8C9+v6fA4nA+0KIy6WU6cAaIAw4S9//Uv37Z1LKeuA/wJOAAD4GfIDXhRDXNPuc2UAf4APgULPXsoGv9McVaFVBL7YSb4MrgUqgGJgGbAdCgHTgFD0uhBDR+jmYBvwB/AZcDCwTQngf5TOUk5BKBMqJ4jb9++NSymuAyYAVOFsI0Q8w66/vBD4CrgWi9G0Nr6WhXdhvA2KBrOYfIoSIAuboT6dJKecDD+jPb9e/v69/v0z/3pAI3hdCCIdY1wJVwC79+S3NPi4dGC+lvFFK+ZPjC1LKNKChCqhYSvk3KeXfmsfbzAop5RzgTf15DdrFviH+kfr3q4BQtPORBRwECoABwJSjfIZyElJtBMqJIk7/vgdASlkohCgEooHeaBfnycCFwFxAAr8KIWYBPwOvoF0AG6pbNgEzgdxWPqdGSpmpP07Wv/fWv38OvARcJIQYAIwBUqWU64QQkUCAvt+1zY7dt9nzDVJK69F+8A7Yo38v1b+nSSntQogK/bm//j1O/z5Q/2orRsUDqBKBcqLI0L8PANDr9yP0bZmAVUp5GRCEdnH7Fe1u+GLAiHaXHoJ2oXsf7eJ9fRuf4+tQZ97f4XOQUpYB3wHBwBv6aw2lhEK0UgDAcCmlkFIKtP+1Mc0+q+4oP7NN/97e/1PbUZ43yNC/f9MQnx5jd7SqM8XDqESguKMnhRDrHL7OABbqrz0ghHgXrV7bBPwipUwFLhdC7EGr378TrQ0AtLvjU4H9aFVGdwOnObzWhJQyH617JsAvQoi3gcf15469dRou/KejlT4+0N8vHWL9WW/A/gStGmhBh84CHNC/9xBCvCWE+EcH39+aj9B+9ll6o/nrQohf9c/r1kmfoZxAVNWQ4o76NXseJqX8VghxKXA/cAlag+jraI3AAClod+PnojUC5wKPAkvQSgF70RqNQ/T9XuPw3Xxz89EuirPQ2gH2Ac9JKT922GcZWgNvNPC7lNKxveEhoAiYh9aAWw5sBT5r588PgJQyQwjxDFrJ5TogCa0R+rhIKXOEEJOAx4BxwES0toKFaOdG8TBCLUyjKIri2VTVkKIoiodTiUBRFMXDqUSgKIri4VQiUBRF8XAnXK+hGTNmyJ9++unoOzooLi4GICwszBkhHTMVV8e4a1zgvrGpuDrGXeOCTolNtPbCCVciKCxUvdsURVE60wmXCBRFUZTOpRKBoiiKh1OJQFEUxcOdcI3FiqJ0DYvFQnZ2NrW1tZ1+bLvdDkBeXl6nH/t4uGtc0P7YfHx86NGjB2azuc39HKlEoChKi7KzswkMDCQuLg5tmYXOY7Vqs2+bTO51CXLXuKB9sUkpKSoqIjs7m/j4+Fb3a05VDSmK0qLa2lrCw8M7PQkoziOEIDw8vMOlOJUIFEVplUoCJ55j+Z15VCKwWu08cMNGvvkgw9WhKIqiuA2PSgRb/yzi87fSee7Bna4ORVGUdjAajYwYMYIhQ4ZwwQUXUFpa6rJYfvvtN9auXdtpx/v222/ZvXt34/OHH36YX3/9tdOO3xEelQhSd5UBUJBbg92u1mFQFHfn6+vLtm3b2LVrF2FhYSxcuPDob3KSthJBQ0NuRzRPBI888ghnnXXWMcd3PDwqEexN0hKB1SrJz61xcTSKonTEhAkTOHjwIAD79u1jxowZjB49mtNPP53k5GRA61o5a9Yshg8fzvDhwxsv3M899xxDhgxhyJAhvPDCCwBkZGQwcOBAbrjhBgYPHsz06dOpqdGuCy+99BKDBg1i2LBhzJ07l4yMDF577TWef/55RowYwerVq5k3bx4333wz48eP5+9//zsLFizgmWeeaYx3yJAhZGRkAPD+++8zbNgwhg8fzlVXXcXatWtZvHgx9913HyNGjGDfvn3MmzePL7/UVkldvnw5I0eOZOjQocyfP5+6Om156759+/Kvf/2LUaNGMXTo0Maf+3i5Xx8pJ2ooEQBkplUSHevnwmgU5cTRV3Rolc12S5OXtWs/m83G8uXLue666wC48cYbee2110hMTGT9+vXceuutrFixgjvuuINJkybxzTffYLPZqKysZPPmzbzzzjusX78eKSXjx49n0qRJhIaGsnfvXj755BPefPNNLr30Ur7++muuuOIK/vvf/7J//368vb0pLS0lJCSEm2++mYCAAO69914AFi1aRHZ2NmvXrsVoNLJgwYIWY09KSuLRRx9l7dq1REREUFxcTFhYGDNnzuT888/nkksuabJ/bW0t8+bNY/ny5fTr14+rr76aV199ldtuuw2AiIgItmzZwiuvvMIzzzzDW2+9dYxn/zCPKRFIKdnrmAj2VbowGkVR2qOmpoYRI0YQHR1NXl4e06ZNo7KykrVr1zJnzhxGjBjBTTfdRG5uLgArVqzglltuAbT2heDgYP744w9mzZqFv78/AQEBXHzxxaxevRqA+Ph4RowYAcDo0aPJzMwEYNiwYVxxxRV8+OGHbfbbnzNnDkajsc2fYcWKFcyZM4eIiAjg6LOHpqSkEB8fT79+2tLd11xzDatWrWp8/eKLL26Mt6HEcbw8pkRQXFBPSVF94/OD+6tcGI2inFjae+feXu2tU29oI6iurubss89m4cKFzJs3j5CQELZt23bccXh7ezc+NhqNjXH98MMPrFq1iu+//57HHnuMnTtb7mDi7+/f+NhkMjWO/gWcMiLbMWbHeI+Xx5QI0pObXvgPZqpEoCgnCj8/P1566SWeffZZ/Pz8iI+P54svvgC00v727dsBOPPMM3n11VcBrTqprKyM008/nW+//Zbq6mqqqqr45ptvOP3001v9LLvdzoEDB5gyZQpPPvkkZWVlVFZWEhgYSEVFRavvi4uLY8uWLQBs2bKF/fv3AzB16lS++OILioqKgMPrCrR2vP79+5ORkUFaWhoAH3zwAZMmTerQ+eooD0oEWlWQr59WjMvJqnZlOIqidNDIkSMZNmwYn3zyCR999BGLFi1i+PDhDB48mO+++w6AF198kZUrVzJ06FBGjx7N7t27GTVqFPPmzWPcuHGMHz+e66+/npEjR7b6OTabjSuvvJKhQ4cycuRI7rjjDkJCQrjgggv45ptvGhuLm5s9ezbFxcUMHjyYl19+ubFqZ/DgwTz44INMmjSJ4cOHc/fddwMwd+5cnn76aUaOHMm+ffsaj+Pj48M777zDnDlzGDp0KAaDgZtvvrkzT+URhJQnVjfKMWPGyE2bNnXoPcXFxTxx926+ee8gQ8aEsmtTCXGJAfyaep6Tomx/XOB+qyGpuDrOXWM7nrj27NnDwIEDOzskwH3n9HHXuKBjsbXyuzt5Vig7Vg1VQwNHhABQmOec+jtFUZQTjfulPSf5zxtDWPl9KSFhXnxtyqCy3EptjRUfX485BYqiKC3ymBJBdA9fRowPx8/fRGi4FwBZqgupoiiK5yQCR2GRWvcrNZZAURTFQxNBaISWCA6kqy6kiqIoHpkIQsK1RHAoW3UhVRRF8dBEoLURHMpWE88pijsLCAhofLx06VL69etHZmbmERO8tWTBggXExsYyYsQIEhMTufjii5vM9nn99dc3ee7JPDMRhGmJIO+gSgSKciJYvnw5d9xxBz/++CO9e/du9/vuuusutm3bxt69e7nsssuYOnUqBQUFALz11lsMGjTouGOz2WzHfQxX88hEEKpXDRWosQSK4vZWrVrFDTfcwJIlS+jTp88xH+eyyy5j+vTpfPzxxwBMnjyZTZs28dprr3Hfffc17vfee+81zvT54YcfMm7cuMbJ7Rou+gEBAdxzzz0MHz6cP//8k0WLFtGvXz/GjRvHDTfc0Pj+goICZs+ezdixYxk7dixr1qwBtNLK/PnzmTx5MgkJCbz00kuNn998yuqG41x66aWccsopTY7TWTyyE31D1VBJQZ2LI1GUE8NjTpqG+sGjTGZXV1fHRRddxG+//caAAQOO+/NGjRp1xBz+s2fPZsKECTz99NMAfPHFFzz00EPs2bOHzz77jDVr1mA2m7n11lv56KOPuPrqq6mqqmL8+PE8++yz5OTkcOWVV7JlyxYCAwOZOnUqw4cPB+DOO+/krrvuYuLEiWRlZXH22WezZ88eAJKTk1m5ciUVFRX079+fW265hdTU1COmrG44zh133MHEiRPJyclpcpzO4JGJICjUC2GAynILFosds9kjC0aK4vbMZjOnnnoqixYt4sUXXzzu47U0pU5kZCQJCQmsW7eO+Ph4UlJSOO2001i4cCGbN29m7NixgDYldlRUFKDN/Dl79mwANmzYwKRJkxqn8ZgzZw6pqakA/Prrr03aIcrLy6ms1Lqtn3feeXh7e+Pt7U1UVBR5eXmtTln966+/kpSUBGiL0zccx7EN5Xh4ZCIwGgXBIV6UFteTk1VN7z6dczIV5WR1tDv3jmrv9MkGg4HPP/+cM888k8cff5wHHnig1X0ffPBBfvjhB4BWp6jeunUrY8aMOWL73Llz+fzzz0lMTOTCCy9ECIGUkmuuuYYnnnjiiP19fHyOug4BaDOZrlu3Dh8fnyNea20K7NaOs2bNGnx8fJwyD5LH3go3VA9l71eDyhTFnfn5+fHDDz80zjjamscee4xt27a1mgS++uorfv75Zy6//PIjXps1axbfffcdn332GZdeeimgTWn95Zdfkp+fD2gT+DUsXONo7Nix/P7775SUlGC1Wvnqq68aX5s+fTr/+9//Gp8fbQ2F1qasnj59epP1mjtjLQZHHp8IDqrpqBXF7YWFhfHTTz/x6KOPsnjxYgAeffRRevTo0fjVkoY1hhMTE/nwww9ZsWIFkZGRR+wXGhrKwIEDycrKYty4cQAMGjSIRx99lOnTpzNs2DCmTZvWuBKao9jYWB544AHGjRvHaaedRlxcHMHBwYC29vGmTZsYNmwYgwYN4rXXXmvz52xtyuqXXnqJzZs3M3LkyHYdp6M8ZhpqgD1b6snPqSEo1ItPXtvHb0tzuf1fg7lzwRBnhNruuE6mqYudyV3jAveNTU1D3THHGldDfb3VamXWrFnMnz+fWbNmuSw2NQ11OzWUCPLU6GJFUY7TggULGDFiBEOGDCE+Pp6LLrrI1SF1iHul4y7UkAhyVSJQFOU4HW2Us7tzaolACDFDCJEihEgTQtzfwuu9hBArhRBbhRA7hBDnOjMeRw3zDeXnqkFlitKaE63qWDm235nTEoEQwggsBM4BBgGXCyGaj+d+CPhcSjkSmAu84qx4mmtYkyA3q5q62hN/iLiidDYfHx+KiopUMjiBSCkpKipqsbtqW5xZNTQOSJNSpgMIIT4FLgQcZ3mSQJD+OBjIOdpBrVZrYwNYe5WWlgJQXWOhzlJLTa2ZwDBJdA9vDmXX8fQDm7j1ocQOHbMzNMTlblRcHeeusR1PXAEBARQXF5OXl9d5AensdjugjRNwJ+4aF7Q/NrPZTGho6BHXybY6DDgzEcQCBxyeZwPjm+2zAPhZCHE74A+c1dKBhBA3AjcCrXYT6yiDQTDnhlj+9690PvxfJhf8JZaeCX6dcmxFORkYjcbGEa6drSFBhYSEOOX4x8pd4wLnxubqxuLLgXellM8KISYAHwghhkgp7Y47SSnfAN4ArfvosXbR8/Otp9Lsha+PVi00ZGQwY88oZ+OqQn7+qoj7nuicJNNR7tblsIGKq+PcNTYVV8e4a1zgnNicWf45CPR0eN5D3+boOuBzACnln4AP4JxbkFYMHaOd1G1/FnXlxyqKorgNZyaCjUCiECJeCOGF1hi8uNk+WcCZAEKIgWiJoMCJMR0hLjEQgOSdZV35sYqiKG7DaYlASmkFbgOWAXvQegclCSEeEULM1He7B7hBCLEd+ASYJ7u4i0Jkdx98/YyUFderMQWKongkp7YRSCmXAkubbXvY4fFu4DRnxnA0BoOgd2IAydvL2Lgqn5l/iXNlOIqiKF3O/fpIuUBcX616aMta1U6gKIrnUYkA6J2orUewc1PHxicoiqKcDFQiAOL0RJCWVK5GUSqK4nFUIgBCI7wJDDZTVWllf2qFq8NRFEXpUioRoK0B2lAq2Li6S3uvKoqiuJxKBLre+niCbetUg7GiKJ5FJQJdnGowVhTFQ6lEoGsYYZyeXIHNZj/K3oqiKCcPlQh0gcFmwqK8qa+zk6Kmm1AUxYO0mQiEEEYhxEwhRP+uCsiV4vpq1UPff5Ll4kgURVG6TpuJQEppAxYBE7omHNcadao28embTyXzzAM7XByNoihK12hP1dBHwDwhxGAhRFjDl7MDc4WxZ0Ry5V/7YjDAa0/sYX9quatDUhRFcbr2JII7gNOBHWhTRBcA+c4MypVOPzuaU6ZEAfDOC6kujkZRFMX52pMIVrXwtdqZQbnaGed0B2DxR5lqYXtFUU56R52GWko5uQvicCtxiQH0iPcne38Viz/OZM78BFeHpCiK4jRHLREIIYKFEO8KIfL0r7eFEMFdEZyrCCGYpJcKFj2boiaiUxTlpNaeqqGXgKuBev1rHvCC80JyD6dMiSQw2Eza7nL++OWQq8NRFEVxmvYkgnOAp6SUPaWUPYGngfOcG5breXkbmXpBDAAvP7LbxdEoiqI4z7GMLPaYepJJ50Tj5W1g85pC9iWrrqSKopyc2pMIlgL3CSGyhBBZwH3AD84Nyz34B5oZPVEbZPbpG/tcHI2iKIpztCcR/A1tUJmv/vUBcJcTY3IrE/QxBd9/koXd7jGFIUVRPMhR5xoC/g94R0oZqX/Nk1KWdkl0biBxSDChEV4UHqrlz5Un7Tg6RVE8WHvmGroI6NMl0bghg0E0jjT+9v39Lo5GURSl8x11QBnwG/CwEMIbyG3YKKX82llBuZuBI0L48YtstXqZoignpfYkgmv17y/p3wVazyGjUyJyQ737BCAEZO2roq7WhrePx/zoiqJ4gPYkgn87PQo35+NnontPP3KyqtmxoYixZ0S5OiRFUZRO02Yi0BuLg4AlUsqVXROSe4pLDCAnq5qNqwtUIlAU5aSiGovbKa6ftqbxtvVqcXtFUU4uqrG4neL1RLB7a4mLI1EURelcqrG4nWJ7+2EyCw5l15CTVUVML39Xh6QoitIp2pMIHsGD5hdqjdFkYPj4cDb/UcjDt27mrSVnuDokRVGUTtGehWkWdEEcJ4RLro1n16Zifvshl2XfHODsWT1dHZKiKMpxa7WxWAixRQgxTQjhry9GM0DfPksI4ZEtpmGR3sy8ojcA9161niTVXqAoykmgrV5DI4BQwAdtMZoYfbsX0K4VyoQQM4QQKUKINCHE/a3sc6kQYrcQIkkI8XG7I3eRqRfEMPq0CGqqbMyb9htFBbWuDklRFOW4HG320WNuG9DHICxEW9hmEHC5EGJQs30SgX8Cp0kpB6PNdOrWDAbB/Lv70WdgICVF9bz0r12uDklRFOW4HK2NYD4wDS0h3CaEuAhIbOexxwFpUsp0ACHEp8CFgONyXzcAC6WUJQBSyqNO72m1Wiku7ljNVGlpKQDVNRbqLLXU1Jo79P6WzJrXjWf+UcFnb6Uz95YYusX6dPgYDXG5GxVXx7lrbCqujnHXuOD4YwsLC2v1taOVCM4GrkfrMnoRcJu+rT1igQMOz7P1bY76Af2EEGuEEOuEEDNaOpAQ4kYhxCYhxKaiIveY+K1HvC8jTw3GapHcdflWdmwsdXVIiqIox6StEsGULvr8RGAy0ANYJYQY2ny9AynlG8AbAGPGjJFtZba2+PnWU2n2wtfH63hibjR7XiJpSdtJS6rkxvM28emqKYw6NbLDxznWn8fZVFwd566xqbg6xl3jAufE1moikFL+fpzHPgg49q/soW9zlA2sl1JagP1CiFS0xLDxOD+7S0RG+/Do66P58JU0Nq4q5N+3b+XbTdMQQrg6NEVRlHY7lsXr22sjkCiEiBdCeAFzgcXN9vkWrTSAECICraoo3YkxdTofPxNX3NqXgCATSVtK+OmrbFeHpCiK0iFOSwRSSitam8IyYA/wuZQySQjxiBBipr7bMqBICLEbWAncJ6V0j0aADvD1M3H+3F4APHL7FspK6l0ckaIoSvu1Z4qJYyalXAosbbbtYYfHErhb/zqhnTEjmj9X5JOZVsntc9Yy9oxITj87mhHjw10dmqIoSptaTQRCiIdbew1ASvlI54dz4jKaDNxwX3/+c+dW1i7PY+3yPL58J53f91/g6tAURVHa1FaJYIHDY4nWhbThMWiT0SkOIrv7ctP9A/n9x1z2bCvlYEa1mqlUURS311YbwRz96xW0QWDXAzfqjxc5P7QT0+BRodz64CD6DdFm4fjj50MujkhRFKVtrSYCKeVXUsqvgPOAF6WUb0sp30Jbl6Arxhic0BIHBwGwdsVRB0sriqK4VHsai72BfwkheqAljmtxbrfTk0LiYK1EsHVtoYsjURRFaVt7EsG9wFtAQ+NxLdocREobevcNwOxl4GBmNUUFtYRHdnwuIkVRlK5w1Dt7KeXHQBzaXEMXAXFSyk+dGtVJwGQ2kNBfW+dYDTJTFMWdtbeKZywwFUgDpgshhjsvpJPHuEnavEMLH0mirtbm4mgURVFadtREIIT4G/A9cDsQDVwMPO3csE4Op57ZjZhefuTn1vL6k3tcHY6iKEqL2lMi+BvwhcPzX4FRTonmJGMwCi69PgGA/z2SxNvPJaMNplYURXEf7UkEocB2h+d+gNE54Zx8Bo4IYeZfeiHt8Pg92/nszRNqTj1FUTxAexLBBuAW/fG9aL2H1jktopPQeXN78Zdb+gDw1P3bqShTk9IpirtI/uoAi0b/TFlmlatDcZn2JILbgRq0KSZmALmcAGsLu5szZkTTZ0Ag5SUWnvzHDurqbLz/UgYfvpyB1Wp3dXiK4rG2vLaPQ1tK2LbIc0vrbY4j0Beg74fWQNxwtUqRUqouMB0khODSGxJ44p7tfPr6PpZ+lkV5qQWA9StLefWbiQSFdM7KaYqitF/ejlIAcjcWM/SI1XQ9Q5slAv2CvwgYI6XcrX+pJHCM4hIDue6efvgHmCgvtRAaYcYvwMj63wp4+T9Jrg5PUTxOdWEd1fl1AKSu8dxZANozsvgjYJ4QYiNatRAAUspip0V1Ehs3KYrBo8PYubGI+EFGstNreeOJDFYvOwTPujo6RfEsBbvKGh8bKqxYqmyY/V3XF0ZKSVVeLQHRvl36ue1pI7gDOB3YARToX2omtePgH2DilCndCAo2E9/fDyEgPbmC2hqrq0NTFI+S75AIBLB7WZ7TP7P8YDV15ZYjtu/7KZdFo3/mxe6LWXbHli7tat6eEsEqDq9BoHQyXz8jsXH+ZO+vYvOaQk47K9rVISmKx2goEdjR7opXv3+Q4RfHdPrn1FdaMPkYKcuq5q3hyzB4Gbjk69MI7uWHb7g3lmorn8/8A7tFa4rd9L+9eAWYOOXe/viGeXd6PM0dNRFIKSc7PQoP12dgENn7q1jzax5hkd70GRiEl5caqqEoztaQCLKBXkDulorOPf7uMlb8Ywf7fswl9pRwQuL8qa/USv4fTl4JQHBvPwb/pTd2i53wgUHEnRXF5pfTWPvEHv58KplJjwzhlPsGsG1hBr3OjCDs1LBOjRHakQiEEAKYCwwFGqbQlFLKezo9Gg/Vd2Agvy/N5d0XUnnjyWTm39OPB54Z6eqwFOWkVV1Ux4HVBY2JYB9aIvCr7ry+MNY6G1/M/IOSfZUAZK8pJHtNIcIAMaeEk7tBa2Yty6xm7RPaFDQJ07sRd2Y3vPxNpH57kKKUCn57cCebFu6lMqeWjB8L6LOmJ9plufO0p41gIVqD8T/Qxg80fCmdpO8gbRGb+jqtWPjlov1qbIGiONG3l//Jl7PWUFdmod6glQgsQJiEvC1lR3t7o4Mbivjk7N/Z9VHGEa9tfDGVkn2V+Hf3YcI/B2Dw0i633ceGcer9A7n4y1OZ8uSwxkWAfcO96DkxAoAep0Yw9anhjLu7HwiozKnFJ9JM7/MiOj0JQPvaCGYBHwOXA3eiTUW9utMj8WBhkT706uNPwaFazF4Gykss/PZjLmdd4Jl9mhXPVZxWQcaKfIZfG4/R7Jz1rw5tK2H/L3kYzAKfEC/SamxYK60ko1V7rH90LwPPim/cP3XxQeorrQy6rCcG4+GYdn2cyZL5G7DV2Un/5RDCKBg8tzcFSWWsfy6F3Z9mATBwTk96TIjA8A8De787yOAregMgDIKwxED6nt+dtO9z6XVGJEbvplXCvSdF4RNipji1kuBRZiIGBjnlnLQnEYSiXfgvR+s++iXwfzRd3F45Tvc/PRxLvZ3li3NY/HEW776QSlFeLZPPjaFbTNd2JVMUV/nxpk1krMinLKOKKY8P69B783eW8t2V6xh0WS9O/efAVu+c1z+bAkDPiZGMu6sfX924CSqt7AIGAgd/LyFjRR5xU7ux6+NMvrtCm1Hnz6eSiRwUREhCANEjQ1h89XqkTRIS70/p/iq+u3I9RrOBn/66haq8WgC6jwklbkoUADFjw4gZe2T9/oj5CfQ4NYKwxIAW4+02PJRuw0MpOFTQofPREe1JBIf0/Q6hrVTmBZQ7LSIPZTQZMJoMjJsUyeKPs1i3Ip91K/JJ6B/It5un4+ffnl+Vopy4akvryfxdu9j9+VQyA+f0IHrkkRdOu9VO8d5KzP5GAmN8MZi0u/Q1j+8hf0cZ+Tt2kvrdQUw+RmLGhdH3vBhMvkYiBgWRt7VEu1M3QN/zuwNgqdeqYWuBPWilgk/PXcXgy3uR9LF2V2/yM5K/vZT87aVNYok7K4rRt/Zl+zv7Sfs+l68uWQtASB9/hl7Vm8ghwQhj21U5wiiI1Je2dZX2XF0eAgqBe4AX0OYdusuJMXm0yO6+jJsUya7NJRgMkJ5SwQPXb+T5j09xSt2gonS2vO0lSDtEjwzt0PvSlx1C2iQIkDbJd1eu5/ptZ7Pro0y2vZVO7sZiQuL9qcqrpVafniWopy/nvjGW7mPDSPk6GwQYvQzk6A2xWasKWPeMVgIwehuQNondKomdEE5oX+0O3KJ32TQaYLMdJo0Ponh9OTvezQCg1xmRjLq5DwfXFVFfZSXr9wJK0iqJGhbMyJv6YDAZGD4/gfKsGvK3l2LyNTL61r6EJQZ2xunsEu3pPvqhw1O1RGUXuO6e/gDkZFXx+N3bWfJpFuMmRfKXm/u6ODJFaVtNcR3vnbYCu9XOrXvPI6inX7vfu3dJDgB9z+1OzoZiCneX896pv5K7qaRxn6IUrXunT5gX9no75Qdq+PScVXQbEYKt3k7EoCDG3dOPgp1lGIwGcjcXU5ZZjd1ip+JgDUiIOzOKETf2abyxsli1YVI+Pkaqqm1EXxJJ4unRlKZXEd4/kOgxoRjNRuLO6gZA4swYStMrCYzxw6TX6RuMggn/GEDSR5lEDAo6oZIAtK/76IoWNksp5ZlOiEdxENPLnyv/2pd3nk/lP3duZeiYMIaO6fw+xIrSWXZ9lImlSusn/9tDO5n53ngK95Tzx3+SGHtnP2LHhwOw7t97Obi6mHG39yc4zp+aonr2/ajNYNPj9AhiTgln1f/t0pKAgAGze5BwdjS1ZRYMRkFIvD9SalNI7/70AHnbSrX3TozAP9IH/6laT/de+nKxoHUZrc6vJbx/EMJwuHTdUCLw1RNBfT30OiOKXme0/DMKIQjtc+SF3ivAxMib+hzfCXSR9lQNTW5hmxpp3EVOmRLF3l1l/PFLHnNPX87/vTiKuTeemH9syslNSsm2tw5P5bzro0yiR4XyxyNJ1BTXU5xWyfwN09jx/n62vpQBwJL5G5scwy/Sm/D+gRiMBhJnxpD2Qy4DL+3JoLla33n/bj6N+wpg0KW9iD0lnM0vp2G3ysaG2Zb4hXvjF950lK7dLrE2lgi0toZaD1xfvD2JINLhcShab6HclndVnGHuTX2oq7OxcVUhD920CZNZcMm1Ca4OS/Eguz7OZPdnWZyxYEiLdf8Fu8vY92Mu+TvKMAeYiBoazME/i/jlb1sb98ndWEzKt9n8dMtmALpPDKYm14pA4BVoQpgEvSdFNnbRHHF9AoMv74XJz9hm+1hwL3+mPjX8mH6uhiRgNAq8zA2JwPPG8LQnETje/ZcDKWiL1dznlIiUI5i9DFx/7wAS+ufw2ZvpLPjrFsZMjCTuBKuHVNzfH4/tJuu3fGZ/dRrrn0th5wcZDLykJ+ueSUHaJft+zGXCPwYy5q99CYj2xVJt5ff/28n651MbrxQx48IYdk0cPiFmqgvr8Y/2pr5ca2T9avYapB0ixwQy9K89CfYPaTMes5N7y9XrPYaMRoHZrCWbujqVCFpSyJFVQSlOiEU5iinndydtTzmb/yjk/us28Okq1UyjHD8pJUhtbv7V/07CbrGzasEuNi9Mw1Zv58+nkgFtTpyyzGrWPLqbdU8lM+a2vqT/nEfBrjKEAaJGhBAY60u/i2LxCfVi1C2HOzcUJZeT9XsB0g7+3bwZfGMshqN0q+wKDe0DJgONJYKG5OBJOjr7qA3IAJ5xVkBK64QQXHFrX3ZuLGbT6kL2p5YT3885Iw0V58jdUsyfTyYz7fmRBHbRQMG8zWUsf2U3kx8bSuQgrb961qp80n7IpSyrmozleRi9DPS7MLZx9ssNz6cCEBDjg7XWTsTAQMbd1Y+ilAqSv8gmb1sp65/T9vHv5sPw6xOIGRfaahVOWP9AQhMDKM+sZtQtffEJcc6o4Y6y6FO5mEwGzPp4hDpVNXQkNfuoe/EPMDFmYgRrl+fz3kt7WfDyaFeHpHTA+mdT2PP5AXzDvTjnlTEdfn/+rlLC+gZi8jGSsSKPjBX5FCSVUZRcjn+UDz3PiKSX/mXyNlJxoIalf9lKbaGF0vRKrtsynbQfcvly1hqkvWlBf/MraQAYzAK7RXtt0GW96O3QABs1NISooSEUJpez890MvEPMjLgh4YhG2OaEEEx+bCiWaiu+Yd5U1LZ/Ph9nslgOtxGoqqE2CCHebuNlKaW8ro33zgBeBIzAW1LK/7ay32y0qSvGSik3HS0mTzfx7GjWLs9n8ceZPPDcCDVl9QmkYbbLvYtzmLFQdmiQ4JbX0/jx5s10HxfG2NsTWXzV+iavFyVXkLWqgDVA99GhXLFiCj9euY3aQm3wVf6OMr6avYb0ZXlIuyRmfBjh/QMJjvdny6v7qM6vwyfMi34XxrDjnQz8o33oMTG8xVgiBgQx5b8dmwLC5GPE5ONef6uNg8mMYNarhurqPa9TZHuqhuahVQ01/MU2f9xiItAXvl8ITEOb3G+jEGKxlHJ3s/0C0SazW3/kUZSWJPQPJKaXHzlZ1Tz/0E7+8dQIAJJ3lOIXYKJXQstzliiuZbfaKUrWBkRVHKyhYFcZUUND2vXeQ1tL+PlOrQdO7oZivr9G+3eJnRBOxOAgQuL8qS2pJ39HGTnri8jdXMJbI5dRml6Fb5SZ/hf2ZNub6aR+pw3a6jExgnF3JWI0axfm0x70ZsPzqfSeEkniBTHY6uyE9Q9sfP1k1aSx2KRd1upViaBFzwDjgEfQpq1+CNjI0UcZjwPSpJTpAEKIT4ELgd3N9vsP8CTt7IVktVopLu7YcsmlpaUAVNdYqLPUUlNr7tD7naW2vvKY33vh1VG89ngGbz6dQvfeRgYMD+KqKdrF4cyZUezcWEZ8f3+e+2QkBoNgw29FvPHkPubdFc/E6ZFtHrvhfLkbd40L2hdbyd4qbA4NkZveTeGUBxOP+r78LWUsvWIbtjo7QX19KU+rQdohpL8fg26PxqjXbfviTej4KCIm+LHxkf2UplchjBA/L4Ruo/1IKI6i6mAtUeODiB4fTLWtUmv1A4zdYcJTWpfkKmsFPS/S2hKcVYVTdRx/+52prKoKAKNRYjBrA+GqauvcpurKUbW1EmOtvcPXvwZhYa0PRm1Pi83VwGdSyhVSyl+Bz4FLpZSbpZSb23hfLHDA4Xm2vq2REGIU0FNK+UNbAQghbhRCbBJCbCoqKmpHyCe//sMCmT1fW1Lvv/cm89Tfk7FZJTar5Oev88g9UMvaX4v46ctcdm0q496rtrFjQxn3z9vB7g7Mt650npJk7eJn8NLuPPcvyW9clzZ3XQm/3rSTQxtLAbBU2yjYXs7ah1P5buYmavLrCe7ny5gH4kiYHYV/jBcDr+vemAQchQ8JoOc07Z8+5rwgwocFIISg32XdGHl3b2JPC23xfZ7IcRyByaxKBG2pAZ4QQpyCViU0Ezjuq7EQwgA8h1b11CYp5RvAGwBjxoyRbWW2tvj51lNp9sLXx+uY3u8svj7HNvPgtJnBZKbWs3FVIdv+LMVkFlz5175k76+irtbG6mV5vPRwGnW1Nmqr7QSFmikvsfC3udt4d9kkhoxu+zwe63l2NneJq7a0nj+fTKa6sI6x/47D5GNsM7akTH0KhVMiyN1cTGlqNclv5GGpsrH2iT1IuyTjxwIGXtKTXR9nNfbgAa0KaMztffEKMBN+TThc03Zs4/8axMCZ1YhIK8IAgcf4N+Zsro7LJLQikdloxM+kjVq2WUwuj6sltaZ6Anz8nfL3355EcD3wIXCV/vyQvu1oDgI9HZ730Lc1CASGAL/pDWbRwGIhxEzVYNx+c2/sQ/KOMipKLYyeGMGEqdrEWDarneTtZRQc0uZFHzwqhPl39+Otp1PYs72Muaev4J1lkxh7etvVRErLCnaX8eHklVQX1AEQNMiLwdf0bPs9SVpJLDjej9gJ4fz5ZDK/PbBTe1FoM2mWH6hhx3sZCAP4R/sQ3NuPhOnRRI8ObTI/ztEIgyC4t79bVnG4k6a9hvRxBHVqiokjSCmXCyF6AwP0TclSyvp2HHsjkCiEiEdLAHOBvzgctwyIaHguhPgNuFclgY4JCDJz6wMD+W1pLhf+pXfjdqPJwFW39+WLt/Yz6rQIpl8ci8lk4LaHB/PuC6lsXF3Io3/bynebp7sw+hOTlJIfb95MdUEdPiFmakstJC3KbjER2Cx2hAEMRgOFeiIIifcnelQYfc/rTtoPuQT18mPQ5T2JGR/Ojnf2U3mwhn6zehA1LFhNPe5kDWsRGA0O3UdVr6GmhBBCauqFEN3R7uC7Ab8c7cBSSqsQ4jZgGVr30bellElCiEeATVLKxZ0QvwIkDAgiYcCRA8v6Dw3hoRdHNtlmMhu48rZEdm4qIWlLCb9+d5Aln2Zx4ZW9mXJeTFeF7DQVuTV8f816xtyeSD8nLfW54939HFhdgFegiUmPDWX5vdspSqqkcGc5YZMOF9vLs6tZNPJnJJL4M7s1TqEcqvfqGnFjAvFndyMwxhej3gV45A1qQsGudLj7qGOJwPPaCFptMRJCLEe/4AshrgOWAk8BPwkhHmrPwaWUS6WU/aSUfaSUj+nbHm4pCUgpJ6vSQNfw8TUyXh8kdMusP1jyaRa3zPqD9b/nt/m+Hz7LYuaoZWTuc48eHy3Zviid/b/ksWT+Riw11mM+Tm1ZPZ9dsJqPzvqN3C3FlKRXsuHFVN4ataxxxsx+F8US1NOvccDVitt3k722sPEY657R2g9qCuvZ/dkB7FZJQIwP3iFaG5UQgpC4gMYkoHS9JlVDDd1HLZ5XImir68AQoKE3z8369/8AvwM3ODMoxfnOODsaACnBy9uA1SK5aeZqcrKqWn3Puy+msntrKa8+3rwHsPtIX3YIgJrCOja+tLfd76vIrWHTwr3UFNdRXVjHx2f9RtqSHDKW5/H26F94pc8P/PK3reRtLcXsZyR+ejf6zdRKUP1mxmDyM1C0s4L3Ji5n37Jcqgpq2famNiXzqJsTGH5dPKNu6cOoe/qzfmMxVqvn3XW6o8YpJgyHZx+1qETQRDBQJIQIBkYCWVLKBcB7QOuTfisnhB7x/ow+LZzQCC/ufmwoQ0aHUllu5f75Gxu7NDqqrbGyS18patlX2dTXu1+DWm1ZPdl/Hu7QtvaJPWx8eS+1pW03adVXWvhk2m8su20Lbw5bxpvDfiJ3Uwl+kd70mhSJwSzwDjYTNSyYETckcM4bYxhzWyJGfXWqgBhfzvhff2Inh4KEJfM38MvftmKpthE5JIiEc7rT78JY+pzTnR//LOLf/9nNyt+ctxC50n6Nk86ZDrcRqEnnmspAW6d4DlrC+Enf3otO6D6quN6N/xiIlNo0B1ffkci/btnM2uV5/Ph5BOde1rS9YPuG4sZ/mooyCyu+z2HG7LZ7yXS1jBX5SJsktE8AJl8DBbvK+fn2Lfz55B4u/vxUekzQ+ibY6m0kf32QlG+yKdxdDlJSkFSOMGgjfkFbfHzsHYmExAcw7u5+AG023HoFmhh8UyyVmfWU7a8i6eMshFGQODOmyfsO5Wm9uNL3VYK+9KHiOk1HFntuiaCtRPB/wAfAULSpqJ/Vt88F1jk5LqWLNFykgkO9uGR+PB+8nMbCR/YyfXZ0k/02rdbuYL19DNTV2vliUbpLEoHdamfPlwfwCjTR5+zuTV5rqBYKHxTIsKvjyFpVwN7FOZRlVvP+xOUMvLQnvuHepH57sPGC38Dka+TUBwZQlFyB3SoZMLtH47w47e25YzAKxt2ZyMp/7sQ33Ith8+LpPrbpIi7l5Vq7RXFxezreKc7W2Ebg0GuoYQ1jT9JqIpBSfqGvV5wA7JFSVgohTGhdQA91VYBK1zn1zG4s+zqb/Jxafvw8l2tua+zdy6Y/tEbQqRfE8NOX2az5JY+yknqCQ7tucF5RUgXLrt1B2b5qEDDnu4mNPYPqyi2kfqsNU+k2MhSjt5H4adH0nhzFjvcySPshl92fHh7oHhDjQ49TI4gYEkR1fh2BPXyJGhJCt+FHrr7VESEJAcz8cDxCgKGF0bsVFdoEcMUlKhG4A8eqIVe3ESxfkYddwrQzu76k2Gb3USllEQ7VQFJKK7Dd2UEprmEwCmbM7sH7/0vjnef2c9WtQzAYBDabnS16b5jxk6NIT6kgZUcZ332YwdW39+uS2KRdsuL2JMr2VeMVYKK+0sq3l/9JzLhwAmN8MZgEVXm1BMf5Ezbo8MptBrOBEdcn0O/CWNJ+yEFKCO8fSPexoU6bUM1obr3prbxCKxGUlVuc8tlKxzQ0FmtTTGi/N6sLSgQ2m+TFl9Ow2yWTTo/Ey6trpwBRE44oTYyfHEVwmIns/TUs/igTgDW/5lFVYSU0wovoHr5MmKr1Ffjo1TSuO3cV91+3ocUG5o6SUpL2Yy65Ww5PqlWWWcXa/+5h9X+SKNxegVeIiekvjyR2QjiWKhuZK/PZ9VFm42jcsOlRzL16I19+ld3k2H6R3gybF8/wa+PpcWqEy2bVbCgRNFQRKa7VcPdvMhkwGsAgwG6ny3t1VVZZsVoldjvk6+1IXcm5C4IqJxyT2cA5l3bj09cO8sS925g2K5aXFiQBMOq0CIQQjJwQwcev7mPfngr27dEGSY06NYJLr0s4rs/e+NLexsXOE8+PYcztiSy9aRNlGYe7tMadH4FvmDfj7+1P9uoCEHDgjyJyNxYTd2Y3chDU19tZtbqAS2b3OK54OpvFYm9cGL2y2ordLjF0YNoIpfM1jCw2mwRCaO0EdfWS2hobAYFdd59c4VBCzMuvo0dPvy77bFCJQGnBuMmhrPm5mAPpNVx+xgp2by3FP9DE2bO0+ngfXyOjTo1g3cp8gsO8KCuu54l7tnHWhbGERbS9UlVLpJTsfD+DX+7SkoDBLNi7JIe9S7S58/27eWOrl/h0M9F7hrZQitFsoLc+r1LvKd2oLa3HHGBi50dZAOQe6vq7qqMpd/hnt9uhtKSesKOs7KU4l2MbAdCYCGpqbAQEdt109Q1VhgAF+vxVXak9K5T1B+4F4tCmigBtZTK1cvpJymAQzJ4fw8sL0tm9tRSACWdGERx2+KI157p4EgYEMnJCOG89k0LKjjL+ce163lh8+hG9bOpqbdTX2QgMPrJhua7Cwjdz/2TfUm1mzr7ndWfAnB7s+fwA+3/OwzfSm9P/NYjAGL82J1Dz0UfrFhRq/0RV1TbKyywEBbvH2hPQ9J8doLCoTiUCF3McWQyHE0JtbdeOk6moPPy3UVjohokA+Bbo32yb5/Wv8jBx/fz49yuj2fRHAWUlFs69tGlX0YAgM5PO0bpvXnVbIo/euZWVS3J569kUrr+nP2l7ykneXsrepHI+eiUNm83OB79OJuvNdFK+ysY/2odBc3uRs6GYfUtzMfsb6XdhLP1nx2I0Gxl1c1+GXRMPSEy+7S+4Ov4T7dtfycgRx9cLqDM1tA80KCioo1/XtLUrrTi8eP3hEgHQWIXXVRyrhtw1EYQBz6PNM6RauDxIRLQPMy45+liByGgfrrq9L28+lcKT923nraeTKcqvwxsYiLZUHcDrp6+guz6hV01xPase3gWAOcDExIcGEjFImwO+vNzCrt3ljB8b1nin1l6FRYe7ZWbsr3KrRNC8RFBU2HldSAsL6/DzM+Lnp2p7O6KlqiGA2i6eeM7xb6PIBWNM2tMa8j7QFwhAKwk0fCkepCKlgvqi1u9UxkyM5IpTwhntbaA4v44eXgZmmwUj0Rai6AF0r7NjA7ymdcPv/O54RXpj8DLgMyWSqgDtAial5PEnk3n08T28+15Gh2KUUja5m8rKqu7wz+lMFc16CnXWP3xRUR033LKZx59M7pTjnUhSUiu4Zv5GtmwtOab3N1QNNYwqbpiBtOY4Jiw8Fo6lRVcMNmzP7cM9aBf+8x22yXa+VzkJVO2vJOnv2/EK82LYwlGYWrjrLNlcjNe6IoYBo3r7Yc2uAZvEu5s3wSNCqKuzkbqumG21dnJ+yWt8n1mA5ftcxJJcpk6OIjExgB07tbaAr749SESENwkJ/gwcEHjEZzZXVmZp0gf8QHZNG3u3rqLCwm+rCjhzSlSn3mFXVGr/7EajwGaTnTaoLGVvJXV1dnbsLMNmkx0uRZ3I1qwtpKCwjiVLchg1suOlv4YpJhpLBCYXVQ05lAhKy7p+jEl7/spXoUoAHq1geT5IqC+qJ3PRfvrcfnjB9dKtJVjLrWR/ktm4zZqp3YkHDQsmano3jPpUDZGnRmBdW4T3vkq8vQxkHqyhptZO9yhv8grrWL4yn+UrtamwY7v5cDCvltff0mbwDA0x07O3N6UlVkpLrfTtE8C/Hx7c5KLXUC1kNgksVknOMfYc+uqbg3z+ZTY7d5bxwP0Dj+kYLWn4Zw8PMZNfVN9pd34NJR+rVZKTU0PPLu566Ep5+VoJMPMYS38NVUOHSwQN6xZ3bWOxY9VQeaW1cQ6wrtKeFcomd0EcipuyW+0U/q7PlCmg4Jc8AgcEEjUtmtxvD5L59v7Gfc2hZrqdE03xn8UEDAwkZGRIkz9mo0EwdWIEUyfqk7/ZJHV1Nnx9TRSX1rP0lzySUiuJjvTmxqt78+emYtIzqiksrqe41EJJ6eE7pS1bS1mxIo9p0w7PidRQLRQb7UNWTg2lZRZq62z4eHds8Fiavt7CH2uL2LK1lFEjQ1rd9/dVBexJKWfuX8KPeife0H00MtyL/KJ6Sks7584v68Dhi+C+tErPSgT64Ku8gjosFntj1U57Ncwr5PLGYoeqIZtNUlZmJSSk63q8taf7qECbaG4o4KNvllLKe5wZmOIeyraWYi2zYA41EzwihMKVBaT/L42Dnx+gLk+78PrE+mKrsRE5NRL/hAD89RW4jsZoFI1VL+GhXlx1aU9KyyyYTQJvbwOTT4tg8mla3X9Wdi2FZeUEBRrJO2RgyS95fPBxFlOmRGHS7+YaSgQhQWZq6uzkFdSxenUB086KbjWGlmRmHr6wLnxlL2+9MabFuzObTfLK6/uoqLAyeIg3I0a2XX3VUCKIivAmKbWS8opOSgQOd8P7M6qY3ClHPTE0lAjsdsjIqCIx8ehViI4Olwiadh+t6eruo806EuTl1XRpImhP+lwIfAT8A/ibw5dykqsvrufABxkA+PcLJPy0CKJmdAMD1OXVIUyC8MmR9JrXm4Rb+xDYwnKZHRUSbMbfv+n9iRCC3j19GTLEn169fZgwNpSIMC8Ki+r59rucxv0aSgQB/kYmjtOWjHz73Qyqq9tu+CsqqqOqStunosJCUXE9ZpMgwM9Ibl4dGRktVzuk7q1o/AfetbPiqD9bQ/G/W6Q2dqC03Np4ITpWNpvkQPbh+DIyW19Y6FiOXdtJVSQWi71TpiFxVFtro8yhPj019ei/g+Yal6rUbyaCg7XSY0Zm13Y0aKw2DNUu/vn5XduFtD2JYBbwsf74TmAl2kplykmqbEMZaQvS2XH7FqozqjEFmQgdpTXEhY4JI+G2vvS8pjcJd/YlYmJEly+wbjQIZkyNBODdDzLYskXrMdIwmCwk0MzoEcH0jPGhrNzKv/+zm9S9LV8k8vJqufHWLVxz3UZ+XZ7XeAGICPOib4I/ABs3Frf43o2bDvdUSUk5+oWjofgfHuJFVIQXNpsk5RguXo4O5dU2mS0z++CxNZC35NXX9/GXq9aTl398o7R3JZUxa85alvyQ20mRaZrHlb6/40mwcYoJvUopsY9W6bFte2mnJ662NJQOo6O0zy/s4tHF7UkEocBq/XEu8CVwo9MiUpxOSkltTg3SfuQfevmuMjKez6JqdxXWCive0d70/EtPvMIPjwo2B5nx6+nXoYFenW3IgCAmnRqO3Q6PPL6HH386RG6udmEIDjFjEIJZ53bHyyzYmVTOXfdub0wYjr765iA1NTaqq2089+Je3ns/A9ATQZyWCLZuK+VgTg37m11oNjkcLyOzFnsL59NRQ4nAP8BInF6Pv31b6TH9/A0aqoV6RGsXkLyCuk6bMO3P9UXU1trZ1EoibK/VfxRit8MPS4+eCIqK6lizthCb7egX4Ty9atKgX8UyMtqfCGw2SUWltTGJeulVQt27e+Hvb6CkzHLE79tZ6uvt1NXZMRi09iPgiI4OVVVWtu+spKrKOVVW7UkEh9DaEg4Bb6EtUKNmLT0BSZukIrmcPQ/tYtvNmznw4eGePrWHajnwYSap/90DNvAf4kuv+XH0mheHV4RPG0d1nbOnRDJ6WDD19Xb+90oaySna3XVokFa8jon24b6/9mXU0CCkhKeeTaFUX7YyL6+W5JQKfv5V68o6aohWrbU7WTtGZJgXffREkJxawZ13b+POu7eRrVfDFJfUk5ZWickoCAowUV8v2ba1kl1JLU+DIaVsLBEE+Jvo3dMXoNX926shEXSP8iYk2IzdDtmtdJstLa2npJ1dVotL6ikp0eLd30rVWGvKyiyN5xlg955yLdbsmiZTKbTkmedTeey/yTz3QupRk0FDiaBXjHYuD3SgNPT2u/u54ur11DWsUKaXCAxCNJYK1qztmoUYG/4ufL2NJPTWbhBW/l7QpAF5564ynvvfQR771z6nxNCeC/pDwD608QS1QBmqjeCEU/hbPpuuWEfS33dQrvfTz1uai91ipzqzil33bOPg5wewllvx7ulF+Hmh+Mb4tri4irswCMGcmTFcemEMEaFmIkLNDO0fSPfow/P3BAaYuOSCGOJ6+lJeYeW5F1LZm1bJDbds5u77tlNfb6dvnB+zzo8hwP9w76Lobj6EhpgJCzVTV2+nutqG1SZ5XV+Q/ht9EZxeMT4kxGn/vE8+mcHf/7mT3bvLqaqykpJa0Vi9kF9Qh90OXmaBl5ehsUSwd1/VcVVBpOt3wRFhXkTrbQ979Auvo5W/53Pt9Zu44ZbNHMxp/YK5aXMJn35+gH3plY3bMo/S7lBQUEdySgXJKRUsWZrLNddt5MZbt1BYVEd1tZX9DnfqmzcVU1xc32LpKftgDdt3lOnxFvDKa2ltfm5DiaBXjC8+3gYqq2ytJkGAykprY9XZH2uLGsecCHG4sRigX6KWWDZuOrIklLq3gm3bS9uMq6MaSoq+PgYSE/zpE+dHdY2N9z84fKO2K0n7nfaIdc5CUO3pPvohgBAiBOgtpez6iTCUY2artZG3NJes9zJAginIhG8vf2pzqrEUW8hdnMOhxQexVljxifUlZGQIXoMFBhfN138sRg0NZtTQ4FZfNxgEcy+K5fnX09m0pZS0fUlYrRJfHwNGg2DKqRGYTYIJY8L4Re8qG6tXtfSN82dDSSl+vkYsFjubt5by8itpLP3pEAYBE8eFUV5tY9uuwxffn385RHmFlXUbipl3dW9iuvvyv1e0i1pkmPaPHBZiJjDAREWllaysanr39u/wz11cUs+69UUIoG+8PyZvA8lplaxYmc855xxexvObb/L59BN9EF89/Oex3bz43Ai8m3Wrraqy8t+nk6mutjFk8OGG/7YSx7ffHeSNRfuP2F5fb+eFF/Yye3YP7A41VR98lMmh51OZOjWMG2+KbfKeZT9rCx/2iPYhN7+WH5flMWN6dxITW+6F1lAiCA0xM6hfIFt2lvHzL4eYf208Foudr785yGmnRdAj1pflK/J4Y9F+qqqsPHT/wCYzfDbv9tsnwQeD0Nocamps+Ppq58lisfPQw7uorrGx6PUxdOvWOSXlhjt/Hx8jQgjOn9aNl97az9KfDhET48tFF8awa7eWIOPinTNJYXu6j8YBXwAjgRlCiP8DfpdSPuyUiJTjIqWkdHMJud8dpPZgLdZyC3a9+Bs6LozIs6IQBkHxn0UULM/ngD6Ng0+sL7GXxmLyN1NHZRufcGIKCTZz3llRfL30EKVlFgIDTNxxXTyBgYf/BU4ZHcLq9UX4ehkau+5NGBNKzqFaTh8XTmFpPb/8XsDSn7QL1qmjQxk4IIjKKisbthYRFGgkeW8tf6wtaux++O77h+/q4nr4MvNsbepsIQRxPX3ZuaeCv927nZHDQ/jbHYls3lKCXcKZU7TFf+x2ycZNxQwaGITNJnnhf3uJj/Pnsjk9Wfx9DhaLpF+8HzExvoSGefHj8nySkitY+Oo+fl2eR3yCD8nJ1QgBk8aHsSO5gqwDNXzwYSbXzY9n67ZSDuXVEh7uTWZmFdXVWtwNd6Cg9W6qrrLi53/k5eLXFdoAwIgwM2aTAbuUDOsfxOqNxWzZXkpxmVZFFBPlTU5+Hbn6Xfzy5cWcNS2M4fr8UhaLvfFYZ50eQfqBalatK+b1N/fx9JPDGjsk2O0SIbTz1zCGICLMi4gIL7bsLOP3VQVcOy+O73/I5b0PM1m3oZgrLu/Fsy/sbYy5YZBia3y8DcR09yE7p5adO0sZN06b+jxpdzmVeh39ipX5XD63F6CVZFJSyjGZDJxxetPOExaLnUN5tfSI9W2yPTOris2bS5g8Kaqxx5CvPvCyezcfZkyN4sfl+bz59n5qa22kpVUiBPSOc041bXta+14DYgEB2NFGGs8FVCJwI1JKqjOqyHgznYpdTasGvLt5EzQihNAxoY1/jEFDgihYoY0Y9u7mTczFWhI4mY0dGcLOPRXs3V/FtIkRTZIAgL+fibtvSgBE411i924+3HZdPAB2KQkNNpO2vxJvs4HpU6MQQqt+uun6aOxS8sJLhyjT/7ED/Y1UVNkwGGDyKeGceUZEYzdFgNHDgklO06aHWLehmGuu20idPtlZaLCZUaNCee+DTL74Kpvhw4Lp2yeADRtL2LCxhB+W5lKvN3ROGB2GEODna2TowCC27izjhx+1htnkZK1+/+xJEUw6LZKBAwJ59d1MfvjxEPUWyRKHBlxDC7WAJpPAapWk7askPiGABY8kMXRIMPOujqOwqI70/VWYTYJbr43Dz6HzQFi4F599l9PY9XbM8BCWrSqgrs5OcKCJsgorr7+WzQ3zfUiI9+f1t9IpK7MQGeZFYt8AevXyY+O2UnYnV/DjT4c495zuFBTU8X8LkggIMPGv/xvEIb1BNSzMi9AQMwH+RgqK6klOruCX5VoJKCW1grff1Uos8T192X+gpnHsQYOWlqaM7+VHdk4tW7cdTgQbHKqKVv9RyOVze7Hy93yeeS6Vhtq9goK6xgWRpJQ88WQy6zYUkxDvz7XXxDFieAivvZnOjz/lYrfDh59kMaC/NvbB1/vwL2DShHD8fI18tSSXDz/JQkqICjfi4+Ocqtr2JIJTgf9yuMvoPrQ5xBQXk1JSvqOMghX5lG4uxqpPambwNRI8PJigIcEY/Y2YAkxHdPE0BZoJnxhBdVY10ed1x+xG8/Y7ixCCqy/rQVGxhahWFtAJamMxEoMQjBoWzKhhLVdDGYRg8IBA1m7UehPNPjeGepsdP28DfRICaN7LdkBiII/8vT9FxRbe/+IA+YX1GAza4KiFr+3jsjk9+EJfcnP7jrLGRtewEDPF+qjkXjE+9O93eBDV+JEhbNXbgM6aGI7By4K/r5GxwyO1O8oefiQm+LM3vYolS3MRAgb1DWBfVjW1dXa6RXhhl1BQpMXSv08ASSkVpKdXkbSnnD3JFexNq2TupT3ZtFn7OXvH+jZJAgAjhwZTVW1jiT6vVGJCAIFBJjKzajjj1HBeeHMfGRm1PKjPQAvg5WXgvDO7aQMNfY2ce2Y3vvohl4Wv7eNQXi0bN5U0jqK+6ZbNVFbZCA40ERJsxmAQjBgSzB/ri3lpYVqTKScyMqsxGQWXnB/DwnczqK7R7uob5nxqSXwvP1avK25SMnLsLpyRVc2OHaW8/Mo+pNTOQebBGt55P4PwMC8mT45k/YZi1m3Qkkf6/ioWPJLEkCHBbN9RhkFoDfy5+XVs2679vnx8mlbVjRkezIYtJRzI0RJedLjzeum158iFwBD9cRRaaSCn9d0VZ7Nb7ZRuLCbn64NUphzuh27wNeAf70/E5Ci8wo7eqBQxKdKZYbols8lAdJTzFoMZPjiItRtLiI70JrGv/1GnnRBCEBHuxV+vjWNHUjm9evjx3ucHyD1Uywv/O9yuUFBcj8Ui6dHdh1vmxVFaZqG0zEJ0lHeTO/nePX05f1o3zAYYOzoMi0Gr5nO8jzzr9Aj2pmsNuGOHBTPr/BgqKi2sWV/CoH4BpKRXseKPQsJDvOgZ60NSSgUpqeVs1RtyrVbJhg3FjY2pfXq3PKXFxPFhhIaYKS+zEBHhRWSkF0MHam0P114Vxdo/y8nLs1GgjwifPaMbA/ofbg8YOzKEymory1YW8OXXWuN8WIiZikpr4wj0i8/p3niOTx0byubtpY1JIKabNzl6VdSAPv6Eh3sxsF8Am/UL7ymjQlizseVZS+N6+iHQLvi/LM/jwIFqDh6swcfbQJ/efiSlVnL/Q1oS6xfvzzVze7Lij0KWry7k6edT+W5JTuO4lknjw7DYJWs3lrB9Rxlms2Du+TEMHhxE8t5KFi87RHGphfBmN2NCCKZMjOD9z7Wbge7dnHez1p5E8CbwmP74I/37/c4JRzma2kO17H5oJ/V68dbgayRocBCBg4PwifXB0FL5XukyvXv4ceNVvQkKMHVoFlBvbyNj9UF7s86N5r3PsgkNNjO4XwCnTwjnmVf3UVNrZ/yIEIxGQXiYF+EtJHshBBPHh7UdY08/Jo4Lo6CwjnPO0qq3ggLNnHOW1i4RGGRi684yBicGNPZEWrWmCMfOTb8uzyNJ72o7qH/r0zoMbuW1qEgzF80Mx5sA7HaJxSLx9j7yb3fKaRFERXiTureC2jo7U0+LoKC0np9W5DNpfDj9+x1OHGEhXlx9aU/e+jATmx3OmRLFDyvyOZRfx1h9XYohAwLZvL2MyHAvzjwjkl0pFfTqfmS9u5+vkW5R3hzKr+P5Fw+3L8T39GPC2DCSUrUEG9PNm4v0ZHTWGRH4+Rr5+bcCUvdqr0eGeTH1jEi8vARmk4HN20s5Z0oUg/XG+AGJAfRN6ENObi0x0UfeoAxIDKBXrC/5hXXE93ROjyEA0Z6ua0KIa4Dz9KdLpJTvOy2ioxgzZozctGlTh95TXKzduezZUk9+Tg1Boc47oR1Roy+96OvTeo8XR9ImSfrnDiqTKzAFmwkcFETYuFBMnby2akNjsTftmzOoq7hrXND5sWn/l6KxOinjQDWZWdVMHB/WpJ3B2XHV19t5/4sDpO3X7rLPGB/GqvWH68pjo7257bqEI6q9nB1XWzKyqsnPr2PMqFAqqywUFtUT18sfg0E7r2s3lBAVrrVFNJ/l0zGuxcsOsXZjCd7eBgb2CaCy2sqUUyPok+BPbZ0dg9Cqs5qrrrGxb38VdXU24nv7N0nYxzKrqMVqx2qF6ooSogcGcuGdY47xzNDqB7er0klK+R7w3rF+unL8bHU2st7eT2VyBUZ/Iz2v6NWu6h/lxNT8YhHX069x7EFX8vIycP0VvSkps1BUVEd8b392pVRQXGrBx9vAhWdHdzgJOFtcLz/iemnnKijQ3KTdRwjBaQ4lprYuypMmhGMyCoYNDKJHrG+T13xaKL008PM1MnRQy/NuHct0LGaTAbMJqo9vNpI2tZoIhBBtjWWWUkq1MI0TSbukdEsJlakV1GRVU5FcgaW4HgRETm1fG4CidJbQYDOheh32KaND+eX3As4/M4peJ/GU18FBZs49q5urw+gSbV3MBdqCNDlAaZdEowBa8THjjXTyms3N4hXuRfjkSALbsVqXojjLGRPCOW1cx9eTVtxXW4ngHWAOEIE26dw7UspfuiQqDyFtEmuFFXOAFVOACVudjbwfcinfVUbpphKEURA4OAjvKG+8o33wjfXF0MGFNxTFGVQSOLm0mgiklNcJIe4ALgXmAz8JIQ4AN0spf2rPwYUQM4AXASPwlpTyv81evxu4HrACBcB8KWXmEQc6CVmrraTcnUpdrtZ1LuH2vlgrrWS9m6HtYICo6d0IGd3xdVgVRVE6os16filllRAiHdgPjEErHbSrXkIIYURb1GYakA1sFEIsllLudthtKzBGSlkthLgFeAq4rOM/hnuTUmItt2Kvt+EdqXVVy3onQ0sCRsAGuYtzMOlD+IOGBRMyOgTf2JO3/lVRFPfRVmPxg8A8IAFYD9wOfCalbG/b9TggTUqZrh/vU+BCoDERSClXOuy/DriyI8GfCKzVVnbdvY1afXRgyJhQfHv6kb/sEBgh6pJwCr4qpiazWmuVMUDElEjMndwlVFEUpTVtlQj+g9ZYnI42ungmMFPv/iSllBce5dixwAGH59nA+Db2vw748WgBW63WxnEB7VVaWgpAdY2FOkstNbVdd5HN+zaf2pxahFmAhNJNJZTqQ9X9x3tjSgTfRB+q99Ro8/709MIeWEcdrpvktZ6uXaavvdw1LnDf2FRcHeOucQFYjTXU2Qwdvv41CAtrfaDh0bqACqCP/uWoU9dwE0JciVb1NKmV129EXxWtR48TZ5ojW62dgh8KAQg7LxTfOB/KVpdhq7PjP8AXcz9tP/+BfloiAHzi3XMRGEVRTl5tJYL44zz2QaCnw/Me+rYmhBBnAQ8Ck1pb60BK+QbwBmgji9vKbG3x862n0uyFr0/n9sGv3FtBRXIF0ed2R+i9Kex1NjI+3o+twoZ3N2/CBkZiMBrwP+fwKOKGUYz+ff0oNBcjrZKQfuF44x7JwB1H8IL7xgXuG5uKq2PcMS6TrQ5vo3+bd/bHfOzWXuiE3jsbgUQhRDxaApgL/MVxByHESOB1YIaUMv84P88lrFVWUh7ZjaXMgsHLQLezo7FWWtl133Zq9dWQQseFYTC23u3T4GWgx196YS2rx9uJE6IpiqK0xGmd0qWUVuA2YBmwB/hcSpkkhHhECDFT3+1pIAD4QgixTQix2FnxOEv2J1lYyrQpgXO+OIC0S/KW5lJ7sAZziJnuF8cS1Mq0xY78evoRNCTkmIagK4qiHA+nThMhpVwKLG227WGHx2c58/OdyVJu4dDiHA4tyQEBBm8Ddfl1FCzP49D32izd4ZMiCWplzhFFURR3oeYLOgY1OTXs+b9d1OvrngaPCMErypuCn/NI1+eQ94r0InCQmgpCURT3pxJBO0ibNgFc4OAgLCUWdv9zB5ZSC15R3kScHk5A/yCkXVJ3qJbynWUgIXhUaJvtAoqiKO5CJYJ2OLQkh8xF+/Ht6Yu0g6XUoi32fkls41oAwiDoPjOGiDMiqM2rJSBRlQYURTkxqERwFNImyV2s1fnXHNB6AXmFexHjkAQcmUO8MIeoKaIVRTlxqETQBmu1ldKNxdQX1GEKNmMKMGEttxA9M0ZNAaEoyklDJYIWSLsk6539WklAH0MdNCyYyEmR2O12tS6woignFZUImrFb7aQ9l0rxH4WNK3yaAk2E6tNBqySgKMrJRiUCB9Im2fe8lgQMXgaiZnQjUB8HYOjAguGKoignEpUIHBz6Poei1VoSiL4whsD+quePoignP3Wb66BwVQEA4WdEqiSgKIrHUIlAV19UR1VaJcIkCBqmpoVQFMVzqESgK9mgLfbg08MXk5+qMVMUxXOoRKBrSAR+cf4ujkRRFKVrqUQA1BXUUra9FIDAgaptQFEUz+LRdSDWaitI2P/qPqRV4pfgj1eYmh5CURTP4rGJwFJmYfutm7FWWAFtPYHIs6LUwjCKongcj60aKlyZ35gEAEInhOMT5R5rBSuKonQljywRSCnJX54HQOS0KPwTA/AKVVVCiqJ4Jo9MBFVpldRkVmPwMRI8IgSjt9HVISmKoriMR1YNFazIB8A/0V8lAUVRPJ7HJQJplxT/WQSgFpZXFEXBAxNBTXolluJ6jAEm/BLU4DFFURSPSwQVm0oA8E/wU4vLK4qi4GGJQEpJeUMi6KtGECuKooAH9RqqPFjLroeTsBTUYfQzEpAY4OqQFEVR3IJHJAK7TbL4ok2UpdcgvAyEnxGJwexRhSFFUZRWecTV0GAUjHuwL8FDgoicFUPomFBXh6QoiuI2PKJEANDnwm7kFQoK0ypcHYqiKIpb8YgSAaAmk1MURWmFxyQCRVEUpWUqESiKong4lQgURVE8nEoEiqIoHk4lAkVRFA/n1EQghJghhEgRQqQJIe5v4XVvIcRn+uvrhRBxzoxHURRFOZLTEoEQwggsBM4BBgGXCyEGNdvtOqBEStkXeB540lnxKIqiKC1z5oCycUCalDIdQAjxKXAhsNthnwuBBfrjL4GXhRBCSilbO6jVaqW4uLhDgZSWlgJQa6+g3lBFdX19h97vLBZRDYBN2l0cSVMqro5z19hUXB3jrnEB1Msq6myGDl//GoSFhbX6mjMTQSxwwOF5NjC+tX2klFYhRBkQDhQ67iSEuBG4EaBHjx7HHFBgNx+KM6qoKbMc8zE6k9Vo1b7b3COeBiqujnPX2FRcHeOucQEYTAb8Qr2dcuwTYooJKeUbwBsAY8aMkW1ltrbMuDYBru3MyI5PQ2Y/1p/HWVRcHeeusam4OsZd4wLnxubMxuKDQE+H5z30bS3uI4QwAcFAkRNjUhRFUZpxZiLYCCQKIeKFEF7AXGBxs30WA9fojy8BVrTVPqAoiqJ0PqdVDel1/rcBywAj8LaUMkkI8QiwSUq5GFgEfCCESAOK0ZKFoiiK0oWc2kYgpVwKLG227WGHx7XAHGfGoCiKorRNjSxWFEXxcCoRKIqieDiVCBRFUTycSgSKoigeTpxovTWFEAVA5jG8NYJmI5bdhIqrY9w1LnDf2FRcHeOuccHxxVYopZzR0gsnXCI4VkKITVLKMa6OozkVV8e4a1zgvrGpuDrGXeMC58WmqoYURVE8nEoEiqIoHs6TEsEbrg6gFSqujnHXuMB9Y1NxdYy7xgVOis1j2ggURVGUlnlSiUBRFEVpgUoEiqIoHu6kTwRCiBlCiBQhRJoQ4n4XxtFTCLFSCLFbCJEkhLhT375ACHFQCLFN/zrXRfFlCCF26jFs0reFCSF+EULs1b+HdnFM/R3OyzYhRLkQ4m+uOGdCiLeFEPlCiF0O21o8P0Lzkv43t0MIMcoFsT0thEjWP/8bIUSIvj1OCFHjcO5e6+K4Wv3dCSH+qZ+zFCHE2V0c12cOMWUIIbbp27vyfLV2jXD+35mU8qT9Qpv+eh+QAHgB24FBLoqlOzBKfxwIpAKD0NZsvtcNzlUGENFs21PA/frj+4EnXfy7PAT0dsU5A84ARgG7jnZ+gHOBHwEBnAKsd0Fs0wGT/vhJh9jiHPdzQVwt/u70/4XtgDcQr//fGrsqrmavPws87ILz1do1wul/Zyd7iWAckCalTJdS1gOfAhe6IhApZa6Ucov+uALYg7Zmszu7EHhPf/wecJHrQuFMYJ+U8lhGlR83KeUqtDUzHLV2fi4E3peadUCIEKJ7V8YmpfxZSmnVn65DWyGwS7VyzlpzIfCplLJOSrkfSEP7/+3SuIQQArgU+MQZn92WNq4RTv87O9kTQSxwwOF5Nm5w8RVCxAEjgfX6ptv0ot3bXV394kACPwshNgshbtS3dZNS5uqPDwHdXBMaoC1a5PjP6Q7nrLXz425/d/PR7hwbxAshtgohfhdCnO6CeFr63bnLOTsdyJNS7nXY1uXnq9k1wul/Zyd7InA7QogA4Cvgb1LKcuBVoA8wAshFK5a6wkQp5SjgHOCvQogzHF+UWlnUJX2NhbbU6UzgC32Tu5yzRq48P20RQjwIWIGP9E25QC8p5UjgbuBjIURQF4bkdr+7Zi6n6Q1Hl5+vFq4RjZz1d3ayJ4KDQE+H5z30bS4hhDCj/YI/klJ+DSClzJNS2qSUduBNnFQcPhop5UH9ez7wjR5HXkNRU/+e74rY0JLTFillnh6jW5wzWj8/bvF3J4SYB5wPXKFfQNCrXor0x5vR6uL7dVVMbfzuXH7OhBAm4GLgs4ZtXX2+WrpG0AV/Zyd7ItgIJAoh4vW7yrnAYlcEotc9LgL2SCmfc9juWKc3C9jV/L1dEJu/ECKw4TFaQ+MutHN1jb7bNcB3XR2brsldmjucM11r52cxcLXeq+MUoMyhaN8lhBAzgL8DM6WU1Q7bI4UQRv1xApAIpHdhXK397hYDc4UQ3kKIeD2uDV0Vl+4sIFlKmd2woSvPV2vXCLri76wrWsNd+YXWsp6KlskfdGEcE9GKdDuAbfrXucAHwE59+2KguwtiS0DrsbEdSGo4T0A4sBzYC/wKhLkgNn+gCAh22Nbl5wwtEeUCFrS62OtaOz9ovTgW6n9zO4ExLogtDa3+uOFv7TV939n673gbsAW4oIvjavV3Bzyon7MU4JyujEvf/i5wc7N9u/J8tXaNcPrfmZpiQlEUxcOd7FVDiqIoylGoRKAoiuLhVCJQFEXxcCoRKIqieDiVCBRFUTycSgSKx9NnmJTNvkqd8DkL9GNf0tnHVpTjYXJ1AIriRraizfQIUO/KQBSlK6kSgaIcVoA2YOdXYLkQYp5+B/+hPhd9oRDi3oadhRA36HPEVwkhNgghJurbvYQQTwghMvW57Fc1+5wpQlsroEAIMUd/z2n6RGy1+vYun/1S8VwqESjKYdPRkkEBTafTmII2Wdoh4GkhxHAhxFS0hcQL0CYj6wUsFkKEo80Zfz/aiNTb0EakOjpTP14w8F9929/RRnj/FXgEKOzsH05RWqOqhhTlsPXAQ/rjEmCo/vhtKeXrQggr8BYwCe3CD/AvKeUvQohewANoC4RcgDZVwGVSm1e+ueeklG8IIW5Bm7sGtOkDzkebUmAL2tQBitIlVIlAUQ4rlFL+qn9tdtgumn13JJt9b4+GRVGsHP4f/AfazJd70ebk2ST05SUVxdlUiUBRDosRQsx1eG7Wv18rhMgC7tCf/442Edg9wL+FEH3QLt4laKuBfQ+MAT4TQnwJDJNS/u0on/1PoA6tOukA2nKNQUDpcf5MinJUKhEoymEjabooyV369+XArUA0cJ+UcjuAvpLb34HngN3AXVLKIiHEfwFf4ApgKu2bTtkO3K5/RhHamrlZx/0TKUo7qNlHFaUV+sIu76Bd/J9xcTiK4jSqjUBRFMXDqRKBoiiKh1MlAkVRFA+nEoGiKIqHU4lAURTFw6lEoCiK4uFUIlAURfFw/w+Fy62CpOM4OwAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Train MSE: 0.2749\n",
      "Test MSE: 0.1889\n"
     ]
    }
   ],
   "source": [
    "import math\n",
    "from tqdm import tqdm\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import torch\n",
    "from torch import nn\n",
    "from torch_geometric_temporal.nn.recurrent import GConvLSTM\n",
    "import torch.nn.functional as F\n",
    "\n",
    "from torch_geometric_temporal.dataset import ChickenpoxDatasetLoader\n",
    "from torch_geometric_temporal.signal import temporal_signal_split\n",
    "\n",
    "loader = ChickenpoxDatasetLoader()\n",
    "\n",
    "lags = 10\n",
    "stride = 1\n",
    "epochs = 200\n",
    "\n",
    "dataset = loader.get_dataset(lags)\n",
    "\n",
    "train_dataset, test_dataset = temporal_signal_split(dataset, train_ratio=0.2)\n",
    "\n",
    "class RecurrentGCVAE_Lin(nn.Module):\n",
    "    def __init__(self, node_features, window_length):\n",
    "        super(RecurrentGCVAE_Lin, self).__init__()\n",
    "        \n",
    "        self.n = window_length\n",
    "        hidden_state_size = 64\n",
    "        latent_dim = 32\n",
    "        \n",
    "        self.encoder = GConvLSTM(node_features, hidden_state_size, 1)\n",
    "        self.latent = nn.Linear(hidden_state_size, 2*latent_dim)\n",
    "\n",
    "        self.decoder = nn.Sequential(\n",
    "            nn.Linear(latent_dim, hidden_state_size),\n",
    "            nn.ReLU(),\n",
    "            nn.Dropout(0.2),\n",
    "            nn.Linear(hidden_state_size, 2*window_length),\n",
    "            nn.Flatten(start_dim=0, end_dim=-2)\n",
    "        )\n",
    "\n",
    "    def forward(self, window, h, c):        \n",
    "        edge_index = window.edge_index\n",
    "        edge_weight = None\n",
    "    \n",
    "        H, C = [], []\n",
    "        \n",
    "        for t in range(self.n):\n",
    "            x = window.x[:,t].unsqueeze(0).T\n",
    "            h, c = self.encoder(x, edge_index, edge_weight, h, c)\n",
    "            \n",
    "            H.append(h.detach())\n",
    "            C.append(c.detach())\n",
    "            \n",
    "        mu_z, logvar_z = self.latent(h).chunk(2, dim=-1)\n",
    "\n",
    "        std_z = torch.exp(0.5 * logvar_z)\n",
    "        eps = torch.randn_like(std_z)\n",
    "\n",
    "        h = std_z * eps + mu_z\n",
    "        \n",
    "        mu_x, logvar_x = self.decoder(h).chunk(2, dim=-1)\n",
    "\n",
    "        kld = torch.mean(-0.5 * torch.mean(1 + logvar_z - (mu_z ** 2) - logvar_z.exp(), dim = 1))\n",
    "\n",
    "        return mu_x, logvar_x, kld, H, C\n",
    "    \n",
    "model = RecurrentGCVAE_Lin(node_features=1, window_length=lags)\n",
    "optimizer = torch.optim.Adam(model.parameters(), lr=0.01)\n",
    "kld_weight = 0.5\n",
    "\n",
    "### TRAIN\n",
    "model.train()\n",
    "\n",
    "loss_history = []\n",
    "kld_history = []\n",
    "\n",
    "for _ in tqdm(range(epochs)):\n",
    "    h, c = None, None\n",
    "    rec_loss = 0\n",
    "    avg_kld = 0\n",
    "    for i, window in enumerate(train_dataset):\n",
    "        optimizer.zero_grad()\n",
    "\n",
    "        mu_x, logvar_x, kld, H, C = model(window, h, c)\n",
    "\n",
    "        log_likelihood = -((window.x - mu_x) ** 2) / (2 * torch.exp(logvar_x)) - logvar_x - math.log(math.sqrt(2 * math.pi))\n",
    "\n",
    "        loss = kld_weight * kld - log_likelihood.mean()\n",
    "\n",
    "        h = H[stride-1]\n",
    "        c = C[stride-1]\n",
    "        \n",
    "        rec_loss += torch.mean((mu_x - window.x)**2).item()\n",
    "        avg_kld += kld.item()\n",
    "                \n",
    "        loss.backward()\n",
    "        optimizer.step()\n",
    "    rec_loss /= i+1\n",
    "    avg_kld /= i+1\n",
    "\n",
    "    kld_history.append(avg_kld)\n",
    "    loss_history.append(rec_loss)\n",
    "    \n",
    "train_loss = rec_loss\n",
    "\n",
    "### TEST \n",
    "model.eval()\n",
    "rec_loss = 0\n",
    "with torch.no_grad():\n",
    "    h, c = None, None\n",
    "    for i, window in enumerate(test_dataset):\n",
    "        mu_x, _, _, H, C = model(window, h, c)\n",
    "        \n",
    "        rec_loss += torch.mean((mu_x - window.x)**2).item()\n",
    "        \n",
    "    rec_loss /= i+1\n",
    "test_loss = rec_loss\n",
    "\n",
    "### RESULTS PLOT\n",
    "colors = ['#2300a8', '#8400a8'] # '#8400a8', '#00A658'\n",
    "plot_dict = {'Reconstruction': (loss_history, colors[0]), 'KL-Divergence': (kld_history, colors[1])}\n",
    "\n",
    "n = len(loss_history)\n",
    "\n",
    "# plot train and val losses and log_likelihood area under the curve\n",
    "fig, ax = plt.subplots()\n",
    "x_axis = list(range(1, n+1))\n",
    "for key, (data, color) in plot_dict.items():\n",
    "    ax.plot(x_axis, data, \n",
    "                label=key, \n",
    "                linewidth=2, \n",
    "                linestyle='-', \n",
    "                alpha=1, \n",
    "                color=color)\n",
    "    ax.fill_between(x_axis, data, \n",
    "                alpha=0.3, \n",
    "                color=color)\n",
    "\n",
    "# figure labels\n",
    "ax.set_title('Loss over time', fontweight='bold')\n",
    "ax.set_xlabel('Epochs', fontweight='bold')\n",
    "ax.set_ylabel('Mean Squared Error', fontweight='bold')\n",
    "ax.legend(loc='upper right')\n",
    "\n",
    "# remove top and right borders\n",
    "ax.spines['top'].set_visible(False)\n",
    "ax.spines['right'].set_visible(False)\n",
    "\n",
    "# adds major gridlines\n",
    "ax.grid(color='grey', linestyle='-', linewidth=0.35, alpha=0.8)\n",
    "\n",
    "# y-axis in log scale (reconstruction loss tends to start high)\n",
    "# ax.set_yscale('log')\n",
    "plt.show()\n",
    "    \n",
    "print(\"Train MSE: {:.4f}\".format(train_loss))\n",
    "print(\"Test MSE: {:.4f}\".format(test_loss))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d27c2cee",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
